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Chapter 1 


INTRODUCTION 


The Geoscience Laser Altimeter System (GLAS) instrument is planned to be launched on 
the Ice, Cloud and land Elevation Satellite (ICES AT) in January 2003 as a part of the 
Earth Observing System (EOS) of NASA. The primary purpose of the GLAS mission is to 
make ice sheet elevation measurements in the polar regions which will be used to determine 
the mass balance of the ice sheets and their contributions to global sea level change [41]. 
In addition, the measurements will meet science objectives to support atmospheric science, 
and land topography application. The laser altimeter will measure the height from the 
spacecraft to Greenland and Antarctic ice sheets to support investigations of the secular 
change in surface elevation, as well as annual, interannual, and other temporal variations. 
To support the science requirements to determine temporal elevation change, the mea- 
surements by the GLAS instrumentation must be very accurate. The ICESAT orbit will 
be near-circular (eccentricity = 0.0013), with a semimajor axis of 6970 km, and it will be 
near polar, with an inclination of 94°, in order to fly over Greenland and Antarctica. The 
orbit is frozen so that the mean perigee is fixed near the north pole. The GLAS/ICESAT 
has a 3 year lifetime requirement with a 5 year goal. The EOS program plans for follow-on 
satellites to extend the science data set to 15 years and longer. 

1.1 GLAS Measurement Requirement 

The GLAS Science Requirements [60] provide the error budget for the instrument and 
ancillary information necessary to meet the science requirements. The most stringent 
requirements are needed for the cryosphere objectives. For example, the accuracy of 
elevation change in the vicinity of the West Antarctic Ross ice streams is 1.5 cm/yr in a 
100 x 100 km 2 area where surface slopes are < 0.6°. In summary, the error budget allows 
10 cm instrument precision, 5 cm radial orbit position, 7.5 cm laser pointing knowledge and 
2 cm or less for other error contributors. The GLAS orbit position will be obtained with 
the Global Positioning System (GPS) and ground-based Satellite Laser Ranging (SLR). 
In addition to the geocentric position vector, the accurate determination of the altimeter 
beam pointing angle in an Earth- fixed Terrestrial Reference Frame (TRF) is also required 
for high-precision satellite laser altimetry. The laser altimeter will provide the range 
between the spacecraft and the illuminated spot on the Earth surface by measuring the 
round-trip travel time of the laser pulse. A one arcsecond error in the laser pointing 
direction produces a 5 cm range measurement error from a 600 km spacecraft altitude on 
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Figure 1.1: Surface slope induced range errors by one arcsecond pointing error 


a 1° slope of the Earth’s surface at the illuminated point(Figure 1.1). Therefore, to be 
able to interpret the laser altimeter height measurements with the required accuracy, it is 
required to determine the laser beam pointing direction to better than 1.5 arcsecond (1 <t), 
which corresponds to a 7.5 cm range error, in post-processing. 

The laser pointing direction is measured with respect to the spacecraft body-fixed 
frame. The GLAS precision attitude determination will first provide the orientation of 
spacecraft body- fixed axes, or attitude, with respect to a Celestial Reference Frame (CRF) 
by star cameras and gyros, to a high degree of accuracy. A specially designed Stellar 
Reference System (SRS) [60] will measure the pointing angle of the GLAS laser beam to 
better than 1.5 arcsecond (la) with respect to the star field for every laser shot fired with 
a frequency of 40 Hz. Calibrated knowledge of the laser pointing direction with respect 
to the spacecraft body-fixed axes, combined with the GLAS geocentric position vector, 
and the measured round trip travel time for the laser pulse enable determination of both 
the spot location on the Earth’s surface and, hence, the surface elevation with respect 
to the adopted TRF. Figure 1.2 illustrates how the GLAS instrument can measure the 
illuminated surface and the surface’s elevation. 

In this document, which is known as the Precision Attitude Determination (PAD) 
Algorithm Theoretical Basis Document (ATBD), we are focusing on the determination of 
the pointing direction from the spacecraft to the laser illuminated surface spot in terms 
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Figure 1.2: Illustration of the GLAS laser altimeter measurement 


of direction cosines (or unit vector). The laser travel time will give the scalar distance 
between the GLAS measurement reference point and the surface spot. Combining this 
scalar distance with the direction cosines gives the laser altitude vector. Since the GLAS 
measurement reference point is not coincident with the spacecraft center of mass from 
which the geocentric position vector is measured, the displacement must be taken into 
account when the GLAS position vector is combined with the laser altitude vector to get 
the laser spot position vector. The detailed procedure for geolocating the spot illuminated 
by the GLAS laser is described in the Laser Footprint Location (Geolocation) and Surface 
Profiles ATBD [42], 

1.2 GLAS Attitude/Pointing Requirement 

Attitude generally refers to the angular orientation of a defined body-fixed coordinate 
system with respect to a separately defined external reference frame, such as a CRF. The 
term, spacecraft attitude, is generally related to the spacecraft body-fixed coordinate frame 
whose origin is the spacecraft center of mass. However, all the instruments related to the 
PAD will be attatched to the GLAS optical bench. Consequently, the more convenient 
choice for the origin of the body-fixed coordinate frame is an instrument reference point 
positioned at the optical bench and the corresponing coordinate frame will define the 
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Figure 1.3: Stellar Reference System on the GLAS 


optical bench attitude as a replacement of the common spacecraft attitude. Thoughout 
this document, we are basically concerned with the optical bench attitude for the GLAS, 
rather than the spacecraft attitude. Spacecraft attitude is sometimes mentioned, but it 
is an alternative terminology for the optical bench attitude when we discuss the GLAS 
PAD. Attitude determination refers to the process of determining the angular orientation 
of the spacecraft-fixed axes or optical bench axes from measurements obtained by various 
attitude sensors. This determination of attitude uses data from appropriate sensors and a 
sophisticated processing of the sensor data. 

Laser pointing , in contrast, refers to the direction of the transmitted laser pulse with 
respect to spacecraft-fixed axes or in an Earth-fixed coordinate system. Pointing determi- 
nation of the GLAS refers to the process of determining the laser pointing at a 40 Hz rate 
in a selected coordinate system. Once the laser pointing is determined in the spacecraft- 
fixed axes, it will be transformed to the vector in the adopted CRF using the determined 
optical bench attitude or an Earth-fixed TRF using the proper transformation between 
CRF and TRF. (Most coordinate systems involved in the GLAS PAD, including CRF and 
TRF, will be described in Section 2.1.) The pointing variation of the laser beam, stemming 
from both the instrument alignment change with respect to the optical bench axes and the 
laser’s own shot to shot fluctuation, will be determined using the SRS specially designed 
for the GLAS PAD (Figure 1.3). 

As will be seen in Chapter 6, the SRS requires the optical bench attitude for the 
determination of the laser pointing vector. In the process, the accuracy of the estimated 
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Table 1.1: Error Budget of the Stellar Reference System 


Error Sources (RSS) 

lcr error 

total Laser Reference System errors 

0.72 arcsecond 

total unmeasured errors 

0.45 arcsecond 

total attitude determination system Errors 

0.85 arcsecond 

TOTAL SRS RSS ERRORS 

1.20 arcsecond 


attitude directly affects the quality of the laser pointing determination. Table 1.1 shows 
the simplified error budget of the SRS [60] when all of the error sources in the system 
are taken into account. Our simulation results, based on the algorithm introduced in this 
document, showed a strong effect of the attitude uncertainty on the overall laser pointing 
error. To meet the laser pointing knowledge requirement of 1.5 arcsecond (lcr), the optical 
bench attitude must be determined with an accuracy of better than one arcsecond. 

Due to the essentially fixed position on the celestial sphere, a star is one of the best 
external sources to provide a reference for attitude determination. The star sensor, the 
attitude sensor which measures relative star positions, currently enables us to determine 
the spacecraft attitude with one arcsecond accuracy level. Among several different star 
sensors, the Charge Coupled Device (CCD) equipped star tracker has been the most recent 
development to provide multiple star images in a single measurement frame. The GLAS 
attitude system includes a Raytheon Optical System Inc.(ROSI) HD-1003 star tracker as 
the primary attitude sensor in order to achieve the required one arcsecond (or better) 
accuracy. Additionally, the Hemispherical Resonator Gyros (HRG) will provide continu- 
ously measured angular rates associated with spacecraft attitude changes. The measured 
angular rates support the attitude propagation between star tracker measurements and 
the attitude prediction when no stars are available. The star tracker and gyros will be 
mounted on a rigid optical bench of the GLAS instrument; thus, this tracker is referred to 
as the optical bench star tracker or the instrument star tracker. 

Two Ball CT-602 star trackers, products of Ball Aerospace, will be mounted on the 
spacecraft structure to support real time attitude determination and control. In the nom- 
inal 600 km altitude orbit, the GLAS Attitude Control System (ACS) must continually 
change the spacecraft orientation with respect to the CRF to maintain a near-nadir laser 
pointing direction. Thus, magnetic torquers and momentum wheels must change the ori- 
entation by 223 arcsec/sec with respect to the stars, which define inertial directions. The 
real time attitude knowledge requirement of the ICES AT is 20 arcseconds (lcr) or better. 
The Ball CT-602 star trackers’ data are intended to be used in the PAD with the optical 
bench star tracker or as a backup unit to the optical bench star tracker. The study in this 
document will only utilize the optical bench star tracker. 

1.3 Outline of Document 

The main purpose of this document is to describe the algorithm that has been developed 
for the GLAS precision attitude and pointing determination, while the title of the docu- 
ment, PAD, only expresses attitude determination. The overview of the process for the 
optical bench attitude determination is illustrated in Figure 1.4 and briefly described here 




based on the figure. At the first step, raw data received from the orbiting GLAS should 
be converted to data written in engineering units. After the conversion, the stars observed 
by the star tracker need to be identified from a star catalog using a priori attitude knowl- 
edge, in order to determine the star coordinates in a CRF defined by the star catalog. 
Stellar aberrations due to the spacecraft velocity with respect to the barycenter of the 
solar system should be corrected before we actually use the identified star coordinates in 
attitude determination algorithms. Modified Fast Optimal Attitude Method (MFOAM) 
and Quarternion Estimator (QUEST) are deterministic approaches for the attitude deter- 
mination using vector measurements. Attitude parameters determined by these methods 
can be used for quality check of the star tracker measurement data. The optical bench 
attitude of ICESat GLAS is determined using star tracker and gyro data, and Extended 
Kalman Filter estimation. 

The next four chapters are related to the attitude data processing shown in Figure 1.4 
except the fact that the conversion of the raw data to engineering units will not be covered 
in this document. In Chapter 2, fundamental coordinate systems, basic attitude equations 
and well-known attitude problems are introduced. The primary attitude sensors for the 
GLAS are the optical bench star tracker and gyro, which are reviewed in Chapter 3. The 
simulation algorithms for star tracker and gyro data are also described in Chapter 3. As 
a part of the attitude determination system, the star catalog is described in Section 3.3. 
The star identification algorithm and simulation results are discussed in Chapter 4. In 
Chapter 5, several deterministic and statistical methods for attitude determination are 
introduced. The simulation results from various methods are presented in Chapter 5. In 
Chapter 6, the implementation of the SRS for the laser pointing determination is discussed 
and some simulation results are presented.* Systematic errors which are imbedded in 
the GLAS PAD system, but intentionally ignored in the previous chapters, will be fully 
discussed. Chapter 7 will address problems related to the actual implementation of the 
attitude/pointing determination system in the GLAS. Quaternions, the preferred set of 
attitude parameters are reviewed in Appendix D as an extension of Chapter 2. 


*An alternative algorithm was used in practice due to degraded performance of the SRS. The new 
algorithm is described in Appendix B. 
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Figure 1.4: Attitude determination overview 









Chapter 2 


FUNDAMENTALS OF 
ATTITUDE DETERMINATION 

2.1 Coordinate Systems 

It is presumed that all of the following coordinate systems have a common origin in the 
spacecraft except for the Terrestrial Reference Frame. The relevant coordinate systems 
are : 


• Optical Bench Coordinate Frame (OBF) 

The OBF is a coordinate frame fixed in the optical bench and is used to define 
the optical bench attitude. The origin of the OBF is an instrument reference point 
located at the optical bench. The orientation of each instrument attached to the 
optical bench will usually be described in terms of the OBF. The OBF x-axis is 
nominally coincident with the optical bench star tracker boresight direction, which 
is zenith pointing. The other two axes complete the proper orthogonal system. The 
adopted OBF for the GLAS is illustrated in Fig 1.3. To maximize the solar arrays’ 
power generation, the ICESAT will be operated in two nominal attitude modes 
[36]. The velocity direction of the ICESAT will change between the OBF ±z-axis 
and the ±y- axis, controlled by yaw maneuvers. Througout this document, zobf 
axis is perpendicular to the orbital plane (upward) and uobf completes the proper 
orthogonal system. Since the definition of the OBF used here is not the same as the 
definition of the planned OBF in Fig 1.3, the appropriate coordinate transformation 
should be applied before actually using the developed algorithm. The angle rotations 
about the (x, y , z)-axes of the OBF are frequently named yaw, roll, and pitch angles, 
respectively. The term, spacecraft attitude, will actually mean the optical bench 
attitude for the GLAS instrument in this document. 

• Celestial Reference Frame (CRF) 

The CRF is a non-rotating coordinate frame defined by appropriate celestial objects. 
Ultimately, it is the inertial reference system which all the other coordinate systems 
are referred to. The simplified description of the CRF can be : the X axis is to 
the vernal equinox direction of a specified date, and Y axis is in the equator and 
the Z axis completes the proper orthogonal system. The CRF will be realized by 
the International Celestial Reference Frame (ICRF) maintained by the International 
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Earth Rotation Service (IERS). 

• Star Tracker Coordinate Frame (SCF) 

The SCF is the coordinate frame fixed to the star tracker mounted on the GLAS 
optical bench. The direction perpendicular to the star tracker field of view (FOV) 
at the center of it is called the boresight direction (BD). While the narrow FOV 
star tracker such as the HD- 1003 can give precise knowledge for the direction of the 
BD vector, it gives relatively poor information about the rotation of the BD vector. 
For the precise determination of the laser altimeter pointing direction, therefore, the 
BD of the optical bench star tracker will be aiming at the zenith direction, which 
is the opposite direction that laser altimeter will be pointing. For our description, 
the zscf is aligned to the BD, the xscf is toward the orbit normal (downward, 
being equal to the nominal —zobf), and the ysCF completes the proper orthogonal 
system (the nominal yoBF )• The orientation of the SCF with respect to the OBF 
will change slowly due to the local deformation of the optical bench and the internal 
error of the star tracker itself. The alignment of the SCF in terms of the OBF 
is assumed to be fixed for the attitude determination process in Chapter 5. The 
alignment variations and the corresponding calibrations will be discussed along with 
the SRS in Chapter 6. 

• Gyro Coordinate Frame (GCF) 

The definition of the GCF is similar to that of the SCF in the sense that the ori- 
entation is defined with respect to the OBF. The GCF may be defined by the axes 
of three orthogonal gyros (usually including a redundant one) in a package, which 
is often called the Inertial Reference Unit (IRU). In this document, we simplify the 
GLAS attitude system by assuming the GCF to be coincident with the OBF. The 
GCF and the OBF may not coincide in the real GLAS attitude system, however, the 
results from this document will not be affected by the change of the GCF orientation 
with respect to the OBF, as long as both GCF and OBF are orthogonal coordinate 
systems. 

• Terrestrial Reference Frame (TRF) 

The TRF is an Earth-fixed coordinate frame whose origin is coincident with the 
center of mass of the Earth. Ultimately, the laser spot location on the Earth’s 
surface will be described in the International Terrestrial Reference Frame (ITRF). 


2.2 Quaternion Representation 

The attitude of the three axis stabilized spacecraft is most conveniently thought of as a 
rotation matrix which transforms a set of reference axes in inertial space to the axes in 
the spacecraft OBF. The rotation matrix is an orthonormal matrix and is also called a 
direction cosine matrix or an attitude matrix. 

Three Euler angles may be used to represent the orientation of a rigid body since the 
rotational motion of the rigid body has three degrees of freedom. There are twelve possible 
sets of Euler angles by the sequence of axes to be rotated. The 3-1-3 Euler angles, which 
rotate the body about the third axis first, the first axis next and the third axis last in newly 
defined coordinate axes obtained by sequential rotations, have been particularly popular 
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Figure 2.1: Coordinate frames 


for attitude determination and control. These Euler angles are shown in Figure 2.1 where 
the coordinate frames are illustrated. 

Although the Euler angles are visually helpful to understand the rotational motion of 
a spacecraft, there is a disadvantage in the Euler angles which is known as a singularity or 
gimbal lock [18]. The 3-1-3 Euler angles are singular when the second Euler angle is 0° or 
180°, because of ambiguous determination of the other two Euler angles. The singularity 
occurs in any sequence of Euler angles and makes Euler angles only infrequently the best 
choice in attitude determination and control application. 

In most modern attitude determination and control, four element quaternions , also 
called Euler parameters , are predominantly used. Their wide use has been the result of 
the following advantages : 

• No geometric singularities 

• Rigorous satisfaction of a set of linear differential equations. 

• No requirement for the evaluation of trigonometric functions 

The lack of trigonometric functions in the computation of quaternions is clearly an ad- 
vantage in time-critical real time operations. Extensive use of trigonometric function in 
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Euler angles will significantly increase the computation time, especially with modest per- 
formance computers such as those used in on-board applications. The quaternions are 
defined based on Euler’s rotation theorem [18] : 

The most general displacement of a rigid body with one point fixed is equivalent 
to a single rotation about some axis through that point. 

For some axis, e, and a single rotation angle, Ad, the quaternions are defined by 


qi 

■ ( A e \ 
= e x sm( — ) 


92 

= e y sm( — ) 


93 

= e z sin(— ) 

(2.1) 

94 

= cost — 
v 2 ' 



where e x , e y and e z are components of rotation axes in terms of the OBF before the 
rotation. Since there are only three degrees of freedom for the rotational motion, the 
following constraint exists in the quaternion representation : 

9i + 92 + 93 + 94 = 1 ( 2 - 2 ) 

The single constraint to be observed is a minor disadvantage associated with the four 
quaternions. The detailed properties of quaternions and relevant equations are summarized 
in Appendix D. 

The quaternion errors 5q are frequently represented by another quaternion rotation, 
which must be composed with the estimated quaternions q in order to obtain the true 
quaternions qtrue as 

qtrue = dq®q, (2.3) 

where the quaternion composition, (g>, is defined in Equation D.9. A benefit of this error 
representation can be seen by applying the small angle approximations to Equation 2.1 : 

T 

x Qy 

4® = y 

= y (2.4) 

5q4 = 1 , 

where 0 X , 0 y , and 9 Z were defined in the previous section as yaw, roll, and pitch angles, 
respectively. Only the vector components of quaternions (see Equation D.l) correspond to 
angle errors and the scalar component becomes insignificant to the first order. By applying 
inverse quaternions (Equation D.2) to Equation 2.3, the error quaternions are expressed 
as 

5q — qtrue ® 9 • (2-5) 
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Even though the Euler angles are not convenient for numerical computations, their 
geometrical significance in illustrating the rotational motion is more apparent than quater- 
nions. Therefore, they are often used for computer input/output and for analysis. In this 
research, simulated attitude data were created by Euler angles and, then, converted to 
the quaternions using Equation D.13. Euler angles can be recalculated from estimated 
quaternions by Equation D.14. 


2.3 Kinematic Equations of Spacecraft Attitude 

If the quaternion composition (Equation D.9) is performed successively in time, the time 
evolution of quaternions during the time interval At is given by 


q(t + At) = q{At) <g> q(t). 


(2.6) 


Let c o x , uj y and ui z be the components of angular velocity vector (cj) at time t, |w| be 
the magnitude of the angular velocity, and A 6 be a rotation angle during At. From the 
definition of quaternions, we can derive 


where 0( cJ) is 


q(t + At) = [ 


P(cD) = 


cos(A0) T | sin(A0) 0(45) 

ri 4x4 T O II 


2 

-*4X4 I 

2 
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0 
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W X 
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Uy 

Wx 

0 

UJ Z 


~Uy 

-UJ Z 

0 _ 


]?(*)» 


(2.7) 


(2.8) 


This equation predicts the attitude at the future time based on the knowledge of the 
current attitude if the axis of rotation is invariant over the time interval At. If the average 
or instantaneous angular velocity of a spacecraft is known during At, the rotation angle is 


A 6 = |w| At (2.9) 

about the rotation axis. Assuming At is small enough, the small angle approximations 

A8 1 , 

sin — « -wAt (2.10) 

lead to the attitude differential equation 


A0 - 
cos — « 1, 


j f q(t) = ifi(w(t))g(t) (2.11) 

from Equation 2.7. Equation 2.11 is the fundamental kinematic equation for the attitude 
determination and can be rearranged in a different order such that [12] 
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Conversely, 


u x = 2(g 4 gi + q 3 q 2 - q-m - qm) 

uj y = 2(-q 3 qi + q A q 2 + q\q 3 - q 2 q±) (2-13) 

u z = 2{q 2 q 3 - qiq 2 + q 4 q 3 - q 3 q 4 ). 

For reference, the 3-1-3 Euler angle representation for angular velocity is 

u x = ip sin 4> sin 8 + 9 cos <f> 

oj y = ip cos 4> sin 6 — 6 siiuj) (2-14) 

(jj z = ip cos 6 + <j), 

where i/j, 6 and <f> are three Euler angles in the sequential order. 


2.4 Dynamical Equations of Spacecraft Attitude 


The rotational motion of a body about its center of mass is 


dA 

dt 



(2.15) 


where T is an applied torque and h is an angular momentum vector. With h described in 
the spacecraft body-fixed axes, it follows that 


hoBF + & x hoBF = 


T, 


(2.16) 


where uj is again an angular velocity vector. Expanding the above equation gives the 
general Euler equations of attitude : 

hx + (dyh z uj z hy = T x 

hy + UJ z hx Ldxhz — ddy (2. IT) 

h z T ^ x dy ujyh x — T^, 

where h x , h y , h z are the angular momentum components along the OBF, and T x ,T y ,T z 
are the body referenced external torque components. For the solution of these equations, 
the external torque of T must be modeled as a function of time as well as a function of 
the position and attitude of the spacecraft. The dominant sources of attitude disturbance 
torques are the Earth’s gravitational and magnetic fields, solar radiation pressure, and 
aerodynamic drag [65]. 

In many spacecraft, gyros are grouped as an IRU. When the gyros are used to measure 
the angular velocity of the spacecraft, the numerical or analytical expression for external 
torques is not necessary. Angular velocities measured by gyros are substituted directly into 
the kinematic equation (Equation 2.11) for attitude prediction. Such gyro measurements 
actually include the effect of all torques acting on the spacecraft. Force model errors 
in this situation will exist only to the extent that the measurements from a gyro unit 
contain errors. Since the HRG will provide the measurement for the angular rate of the 
ICESAT, the Euler equations are not required for the attitude determination. The attitude 
determination in the event of gyro failure will be discussed in Chapter 5. 
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2.5 Attitude Determination Problem 


The minimum required knowledge for three-axis attitude determination is the direction 
vectors to two celestial objects which are represented in the OBF (or the Spacecraft Body- 
Fixed Coordinate Frame generally) and are also known in the reference frame, such as the 
CRF. Since the stars are measured in the SCF, the unit vectors in the OBF are determined 
using the rotation matrix between two coordinate systems. Denote the two unit vector 
measurements by W\ and W 2 in the OBF and V\ and V 2 in the CRF. To obtain V\ and 
V 2 , the measured stars in the SCF must be identified in a given star catalog with the 
star identification algorithms developed in Chapter 4. A simple attitude determination 
problem is given as : 


Find an attitude matrix A, to satisfy 

W] = AV i, W 2 = AV 2 , (2.18) 

where the measured vectors and the catalog vectors require 

W 1 -W 2 = Vi-V 2 (2.19) 

within the measurement and the catalog error bound. 

A simple algorithm to find the attitude matrix from any two vector measurements 
is called TRIAD [49] or an algebraic method [65]. The method has been applied for 
at least three decades and was employed in several missions, usually for coarse attitude 
determination. Whereas it is relatively easy to evaluate the TRIAD attitude matrix, it has 
many disadvantages. The most serious disadvantage might occur when more than three 
unit vectors are observed simultaneously, which is a common case in CCD star tracker 
measurements. (In Section 4.1, three observation vectors are the minimum number for 
utilizing a pattern matching algorithm.) 

To take advantage of multiple unit vectors simultaneously obtained by a CCD star 
tracker (or the combination of multiple sensors), a least squares attitude problem was 
suggested in the early 1960’s by Wahba [64], The well-known Wahba Problem is : 


Find the proper orthogonal matrix A that minimizes the loss function, J(A), 

1 n 

J (A) = -Y J a i\W i -AV i \ 2 , (2.20) 

i = 1 

where the unit vectors V are representations in a reference frame of the di- 
rections to some observed objects, the W % are the unit vector representations 
of corresponding observations in the spacecraft body frame, the ai are positive 
weights, and n is the number of observations 

The Wahba Problem is basically a weighted least squares problem for the attitude matrix, 
A. It is also known to be equivalent to a maximum likelihood estimation problem for a 
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simple, but realistic probabilistic model for vector measurements [46]. For error-free ob- 
servations (and also error- free catalog positions), the true attitude matrix At rU e will drive 
the loss function, J(A), to be zero. In a practical situation, the A must be found that min- 
imizes J(A ). The solutions of the Wahba Problem, which are deterministic approaches to 
the computation of the attitude matrix (or quaternions) will be introduced in Chapter 5. 
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Chapter 3 


MEASUREMENT SYSTEM 


Many spacecraft use gyro units to continuously measure their angular velocities. Tradi- 
tional mechanical gyros react to the motion of the host spacecraft based on the principle of 
conservation of angular momentum. Non-mechanical gyros have been constructed on phys- 
ical phenomena such as general relativity or the inertial vibration property of a standing 
vibration wave on a hemispherical body. Such gyros are usually sensitive to high frequency 
noise and able to measure attitude change very accurately. However, slowly drifting gyro 
biases will produce deviations of predicted attitude from true attitude. Some external 
sources such as the Sun, the Moon, the Earth and stars must be observed in order to 
prevent gyro biases from deteriorating the attitude determination process based on gyros. 
Measurements from external sources are relatively insensitive to the high frequency of at- 
titude change due to instrument noise and jitter. However, such measurements provide 
good information in the low frequency of attitude motion because the positions of celestial 
objects are well-predicted (Sun, Moon, and Earth) or they are essentially fixed in space 
(stars). Therefore, gyros and sensors for the external sources are generally used together 
in the attitude determination system. 

The CCD star tracker which measures multiple stars with a frequency of 10 Hz also 
enables angular rate information to be inferred and may be used alone for the accurate 
attitude determination. However, in reality, the high frequency jitter in the spacecraft 
motion, the irregular distribution of stars, CCD limitations for data read-out and the 
blockage of star measurements due to the Sun or the Moon require measurements of the 
spacecraft angular velocity from gyros for continuous attitude determination. 

This chapter introduces CCD star trackers and gyros which will be used in the ICE- 
SAT/GLAS. In addition, a star catalog which is an essential component of the attitude 
determination using star sensors will be discussed. 

3.1 CCD Star Tracker 

The CCD image detector was developed in 1970 by Boyle and Smith at Bell Laboratories 
[10]. Unlike the traditional light detectors, two-dimensional CCD detectors allow one to 
obtain images of several objects in a single exposure and provide much better photometric 
accuracy than photometric emulsions do. In structure, the CCD is a two dimensional 
array of photo-detectors, and each individual detector is described as a picture element 
or pixel (see Figure 3.1). After a period of being exposed to light, called the integration 
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. star 



Figure 3.1: Simple illustration of the CCD star tracker 


time, the photocharges in each pixel are transferred to an output stage by the external 
manipulation of voltages. The charges are measured, digitized and read into a computer’s 
memory one pixel at a time, row by row. The area viewed by a two dimensional CCD array 
can be quickly reconstructed as a digitized image in a computer in order to analyze or 
process the light intensity distribution of an extended field. The high quantum efficiency 
of the CCD allows it to record about 60% of the photons falling on each pixel. The 
reconstructed image is extremely similar to the area projected onto the CCD due to the 
direct relationship between the exposure and the intensity of the recorded image over 
a broad range of exposures from threshold to saturation. Since each pixel represents a 
specific location on the CCD, a computer can be programmed to perform an analysis of a 
star field automatically. More detailed knowledge of the CCD can be obtained from many 
sources [10] [21] [57]. 

In space applications, CCDs have been used on many imaging missions including the 
Hubble Space Telescope. The CCD star tracker was developed in the late 1970’s and it 
has recently begun to be used in space missions as a component in state-of-the-art attitude 
determination hardware. The traditional star sensor detects only one or two stars in its 
FOV, and therefore, has been used with other types of attitude hardware like Sun sensors, 
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horizon sensors or magnetometers. Since the CCD star tracker observes multiple stars 
simultaneously, it is sometimes referred to as a star camera and some spacecraft use a 
CCD star tracker as the sole attitude sensor except for initial attitude acquisition in real 
time applications. Even though the initial attitude acquisition may require other sensors 
for coarse attitude determination, this may be also performed by a star tracker with a 
wide FOV [33] [24], A simplified illustration of the CCD array and lens in the star tracker 
is shown in Figure 3.1. The following sections describe characteristics of commercially 
available CCD star trackers. 

3.1.1 HD-1003 Star Tracker 

For the PAD of the GLAS, the Raytheon Optical System Inc.(ROSI) HD-1003 [11], a 
CCD star tracker, will be mounted on the instrument optical bench. It is an electro- 
optical sensor that implements CCD array detectors to search for and track up to six stars 
in an 8° x 8° FOV with an array of 512 x 512 pixels. The star tracker operates at a 10 
Hz rate, thereby measuring coordinates of tracking stars as well as their light intensities 
every 0.1 second. It provides the position of a star with a six arcseconds (la) error in 
each axis of pitch and roll. A magnitude measurement is given within 0.2 magnitude 
(ler) uncertainty. The nominal sensitivity range of the star magnitude is between 2.0 and 
6.0. Functionally, the HD-1003 star tracker is operated in search and track mode. The 
star tracker searches the entire FOV to find the six brightest stars in search mode. After 
acquisition, it continues tracking these stars and periodically computes updated angular 
positions as the stars pass across the FOV in track mode. The performance characteristics 
of the ROSI HD-1003 are summarized and compared to the characteristics of the Ball 
CT-601 and the Lockheed AST in Table 3.1. 

The HD- 1003 star tracker manufactured by Raytheon for GLAS was tested during 
the period between summer 1999 and spring 2000. Various tests were conducted, such 
as mechanical properties measurement, spectral calibration, thermal vacuum segment and 
acceptance vibration. Table 3.2 shows the results of the final performance test which 
was performed in an air conditioned room after electronics and software upgrades [13]. 
The static accuracy test measures the position accuracy of stars that are fixed relative 
to the tracker. The dynamic accuracy test measures the star location accuracy while 
tracking a moving star at a rate of 0.20 deg/second. The table presents only the post- 
vibration results, but the pre-vibration results showed similar values. All requirements 
were reported to have been met based on the pass/fail criteria specified in both pre/post 
vibration. The overall performance of the HD- 1003 was not compromised by exposure 
to the random vibration, thermal or vacuum environment. After the final performance 
test, point sources were generated by the Scene Simulator computer in order to simulate 
multiple stars for the tracking performance test. For the Scene Simulator test, six stars 
were simultaneously placed in the full FOV (8 x 8°) and one star in the reduced FOV 
(0.5 x 0.5°). No test failures or anomalies were encountered during the testing. 

3.1.2 Star Measurement on CCD Star Tracker 

In the focused image on a CCD array, a star will appear as a point source because all 
luminous power from the star will end up in one pixel. By slightly defocusing the image, 
however, a star will cover several pixels, usually in an area of 3 x 3 or 4 x 4 pixels [59]. 
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Table 3.1: Characteristics of the commercial star trackers [11] [7] [62], The performance 
parameters of the Ball CT-602 is identical to the CT-601 except the weight and the power 
requirement. 


CHARACTERISTICS 

HD-1003 

CT-601 

AST f 

Field of View 

00 

o 

X 

00 

o 

00 

o 

X 

00 

o 

8.8° x 8.8° 

Sensitivity Range (M v ) 

+2 to +6 

+1 to +6 

+7.5 

Accuracy (arcsec, la) $ 
in Roll and Pitch 

6 

3 (bias) 
5 (random) 

1 

Update Rate (Hz) 

10 

10 

10 or 5 

Acquisition Time 

(Full Field, second) 

6 

5 

N/A 

Maximum Number 

of Stars Tracked 

6 

5 

20 

Maximum Power (W) 

12 

(average) 

<10 
@28V d.c. 

4(0°)/7(50°) 

Maximum Weight 

8 lb 

(with shade) 

18 lb 
(w/o shade) 

3.5 kg 
(w/o shade) 

Operating 

Temperature (°C) 

-20 to +60 

-30 to +50 

-30 to +50 


f development goals 
J Yaw is worse 

The size of the illuminated square area will be determined by the brightness of each star. 
From the illuminated pixels, the centroid is mathematically computed with an accuracy of 
0.1 pixel level [21]. The summation of transferred charges on those pixels is proportional 
to the brightness of the star, from which the instrumental magnitude of the star can be 
computed. Only the pixels that receive more than a certain amount (threshold) of photons 
will be recorded. Dimmer stars will be difficult to detect due to the increased noise. If 
the star is too bright, saturation may occur due to overflow of the photons to adjacent 
pixels. The sensitivity range of star magnitude is between 2.0 and 6.0 visual magnitude 
for the HD-1003. The central position of the star image can be determined precisely using 
the techniques such as Image Moment Analysis (up to the level of j^) or 1-Dimensional 
Marginal Fitting (^ pixel or better). The Point Spread Function [21] is the method to 
discriminate the overlapped stars in some pixels. 

The two angles of star position measured by the star tracker are converted to a unit 
vector, SscFi in the SCF, where the BD is regarded as the third axis. Then, the coordi- 
nates in the SCF are transformed to those in the OBF using the transformation matrix, 
M : 

SsCF = MSoBF, (3.1) 

where Sqbf is the unit vector of the same star in the OBF. The transformation matrix M 
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Table 3.2: Test results of HD-1003 star tracker from final performance tests. The specific 
requirements were given for GLAS. For the dynamic accuracy test, the star position error 
is effectively rms of bias and random noise. 


Star Position Error 


Error 

Type 

Source 

Magnitude 

Requirement 
(la per axis) 

Measured ROW 
axis error (la) 

Measured COL 
axis error (la) 

static 

3.0 rrii 

< 4.5 arcsec 

1.23 arcsec 

1.28 arcsec 

bias 

5.5 mi 

< 4.5 arcsec 

2.44 arcsec 

2.60 arcsec 

static 

3.0 rrii 

<3.3 arcsec 

0.94 arcsec 

0.54 arcsec 

random noise 

5.5 mi 

<3.3 arcsec 

2.53 arcsec 

2.16 arcsec 

dynamic rms 

5.5 mi 

<5.0 arcsec 

2.95 arcsec 

3.23 arcsec 


Magnitude Error 


Error 

Type 

Source 

Movement 

Source 

Magnitude 

Required 

Measured 

error 

30 - 

static 

3.0 mt 

< 0.24 mi 

0.12 mi 

repeatability 

static 

5.5 mt 

< 0.24 mi 

0.12 mi 


dynamic 

5.5 mt 

< 0.24 mi 

0.18 mi 


static 

3.0 mt 

± 0.12 mi 

-0.02 rrii 

bias 

static 

5.5 mi 

± 0.12 mi 

0.01 mi 


dynamic 

5.5 mi 

± 0.12 mi 

0.05 mi 


Boresight (Pitch & Roll) and Yaw Measurement Results 


Measured 

Required 

Measured 

Required 

Measured 

Axis 

Range 

Position 

Knowledge 

Knowledge 

Pitch 

255.5 ± 3 pixels 

254.780 

± 2 arcsec 

±1.6 arcsec 

Roll 

255.5 ± 3 pixels 

255.385 

± 2 arcsec 

±1.6 arcsec 

Yaw 

0° ±0.20° 

-0.004° 

±10 arcsec 

± 3 arcsec 
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is determined in prelaunch instrument calibration, but it may be changed slightly by the 
launch loads and other effects like temperature variation and atmospheric drag in orbit. 
This is known as the star tracker alignment error that is also represented by the star 
tracker BD excursion. 

Using the SRS, the GLAS laser beam pointing direction would be determined against 
the star field viewed by the star tracker that defines the SCF. In other words, the laser beam 
pointing direction can be directly associated with the SCF through the SRS. Therefore, 
knowledge of the orientation of the SCF, not the orientation of the OBF, in terms of the 
CRF is needed for the GLAS PAD. This fact will allow us to determine the orientation of 
the OBF with respect to the CRF under the assumption that the star tracker is rigid in 
terms of the OBF. This assumption will not restrict the implementation of the estimated 
attitude to the laser pointing determination, as will be described in Chapter 6. 

3.1.3 Star Tracker Data Simulation 

The orbit of the ICESAT is near-circular (eccentricity = 0.0013), with a semimajor axis 
of 6970 km and an inclination of 94°. The BD is assumed to look at the zenith which is 
opposite to the laser pointing direction. The misalignment of the BD to the laser pointing 
direction will be discussed in Chapter 6. The BD can be calculated in the CRF along 
the ICESAT orbit at a specified interval, such as 0.1 second. The 8° x 8° FOV whose 
center is the BD is constructed by specifying the FOV boundary in right ascension (a) 
and declination (5). The stars located in a FOV are found in a star catalog and only the 
brightest stars (up to six) are selected as measured stars by the CCD star tracker at each 
measurement time. The star positions (in a and 5) and star magnitudes are given as the 
input parameters for the simulation. A schematic of the procedure is given in Figure 3.2. 

The stars are selected from the star catalog in which star positions are described in a 
CRF. Since the real star tracker measures stars in the SCF, it is necessary to transform 
the coordinates of the selected stars. The relations between coordinate frames, described 
in Section 2.1, are used for coordinate transformations. In our simulation configuration, 
it is assumed that the OBF is obtained by the 3-1-3 Euler angle rotation from the CRF 
(Figure 2.1). The Euler angles are the longitude of the ascending node, D, inclination, i , 
and the argument of latitude, u. These 3-1-3 Euler angle rotations align xobf with the 
BD, zobf with the orbit normal direction (upward), and yoBF with the velocity vector 
direction. The SCF is obtained by simple rotation about yoBF to make the zscf to be 
the BD of star tracker. 

Thus, the selected stars from the star catalog, whose coordinates are given in terms 
of the CRF, are now transformed to the OBF and the SCF successively by the following 
procedure : 

1. a and 5 to unit vector 

First, the two angles in the CRF are converted to a unit vector in the CRF by 

X cos 8 cos a 

Y = cos 8 sin a . (3-2) 

Z sin 8 

2. CRF to OBF 

The unit vector of the star in the CRF is transformed to the unit vector in the OBF 


24 



Start 



Set Boresight, Direction 



Figure 3.2: Flowchart of the star selection process for star tracker measurement data 
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through the 3-1-3 Euler angle rotations : 


XOBF 

UOBF 

ZQBF 


cos u 

sin u 

O' 

— sin u 

cos u 

0 

0 

0 

1_ 

cos n 

sin Q 

0 

— sin 12 

cos 12 

0 

0 

0 

1 


1 0 0 
0 cos i sin i 

0 — sin i cos i 

X 

Y . 

Z 


(3.3) 


The magnitude of 17 was varied in the simulation to comply with the 183 day repeat 
ground track of the ICESAT. 


3. OBF to SCF 

The unit vector in the OBF is transformed to the unit vector in the SCF by the 
90° rotation about yoBF- This corresponds to the rotational matrix M defined in 
Equation 3.1. Then 


'xscf 


1 

H 

1 

o 

o 


XOBF 

yscF 

= 

0 1 0 


yoBF 

_ZSCF_ 


O 

O 

i— l 


ZQBF _ 


This rotation aligns zscf with the BD, xscf with orbit normal direction (downward) 
and ysCF completes the right hand coordinate system. 


4. Unit vector to position angles (4> and A) 

The two measurement angles, (j) and A, determined by the CCD star tracker are 
computed from the unit vector in the SCF. First, the focal length / (Figure 3.1) of 
the star tracker optical system is given in millimeter (mm) unit. If x mm and y mm 
are the distances to the star image from the center of the CCD array using the same 
mm units as /, the star position on the CCD detector array is given by 


%mm 


xscf/ zscf X / 

Vmm 


yscF/ ZSCF X / 


Then, the two measurement angles <f> and A in the SCF, analogous to a and 5 in the 
CRF, are defined by 


%mm 


Vmm 


tan (j) X / 
tan A x 


/ 


cos(4 > ) 


(3.6) 

(3.7) 


From the known noise characteristics of the CCD star tracker, noise components are added 
to 4> and A. A measurement noise is also added to the magnitude. A set of realistic noise 
characteristics for the GLAS star tracker are given in Table 3.3. The noisy position angles 
and magnitude in the SCF are considered to be the simulated measurements obtained by 
a star tracker. The purpose of the attitude determination that will be discussed in later 
chapters is to find input Euler angles which were used for data simulation. As mentioned 
in Chapter 2, quaternions will be used instead of Euler angles. The flowchart showing the 
procedure to get simulated star tracker data from the stars which were selected in a star 
catalog is illustrated in Figure 3.3. 
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Figure 3.3: Flowchart of the star tracker data generation from the selected stars 
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3.2 Hemispherical Resonator Gyro (HRG) 

Gyros provide the spacecraft angular rate, even though several independent technologies 
are used to design different gyros, such as the traditional mechanical gyro, laser gyros 
and the HRGs. The ICESAT will use a set of HRGs mounted on the instrument optical 
bench. An HRG is based on the inertial sensing property of a standing vibration wave on 
a hemispherical body [66]. It was noted that the location of a vibration pattern at the 
rim of a hemispherical shell precesses relative to a reference on the shell when the shell 
is rotated about its axis of symmetry. The vibration pattern precession was observed to 
be a constant fraction of the inertial input. Measuring the amount of precession provides 
the inertial rotation of the input axis. A carefully designed IRU consisting of a few HRGs 
provides the inherent small size, high reliability characteristics and precision performance 
with relatively long life time expected in space applications. 

3.2.1 Gyro Model 

Farrenkopf’s gyro model [14] [65] has been widely accepted and used for many years. It 
basically separates the gyro noise into three noise types called electronic noise, float torque 
noise and float torque derivative noise. The electronic noise is a kind of scale factor error 
generated from electronic part of gyro. It is a colored noise, but it can be treated as white 
noise if the gyro time constant is much smaller than the gyro readout time interval, which 
is usually true for the modern gyros. The float torque noise (rate white noise) is simple 
white Gaussian noise superimposed on the gyro drift rate. The float torque derivative 
noise (drift rate ramp) is the cause of the gyro drift rate bias since its integration appears 
to be rate random walk. This gyro characteristic eventually causes the measured data to 
deviate systematically from the true angular rate and this is why gyros need help from the 
external sources such as star, Sun or Earth. The gyro rate bias, due to the float torque 
derivative, as well as the attitude parameters (e.g., quaternions) can be estimated with 
appropriate estimation algorithms. 

3.2.2 Gyro Data Simulation 

To keep the laser altimeter pointing direction to the nadir direction, the ICESAT must 
rotate about the orbit normal axis ( zobf ) once per orbital period. The nominal (and 
unperturbed) angular velocity with respect to the OBF can be assumed to be 

uj x 0.0 

UJy = 0.0 , 

J T - 

where T is the orbital period of the ICESAT. However, when the attitude determination 
process is simulated for the multiple orbits of the ICESAT, the simple angular velocity in 
Equation 3.8 makes the CCD star tracker to repeatedly image the same group of stars. To 
avoid this unrealistic situation, Equation 2.14 is used to provide nominal angular velocity, 


(3.8) 
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based on 3-1-3 Euler angle rotations from the CRF to the OBF : 

oj x = Q 

UJy = Cl 

uj z = Cl 

The expected value of Cl is approximately 0.5 ° /day for the ICESAT [26] and is zero by 
assuming the inclination of the spacecraft to be constant in the simulation. 

To generate noisy gyro data, the expected gyro noises must be added. Assuming 
that the direction of the gyro’s input axis is aligned with the OBF, the simulated (gyro- 
measured) angular velocity of the zth gyro, gi , is obtained from Farrenkopf’s gyro model 
[65] : 

9i = (1 + hi) + bio + bi + 77 * 1 , (3.10) 

where kj is a scale factor error, cjj is a true or nominal angular velocity, 6 ,;o is an initial bias 
error, bi is a gyro bias which is slowly varying in the orbit and rjn is a white noise on gyro 
rate. The scale factor relates the gyro output counts to the physical unit measurements and 
is a function of the angular rate. Assuming that the kiUi term is negligible in comparison 
with 6 *o and bi, Equation 3.10 reduces the gyro noise into white noise and random walk 
(non-white noise) components. 

If the gyro measurement vector, g, is measured in the GCF while other vectors are 
described in the OBF, then the vector form of Equation 3.10 becomes 

g = G (u) + bo + b + rfi), (3-11) 

where u is the nominal spacecraft angular velocity vector in the OBF and G is an orthog- 
onal matrix transforming the OBF into the GCF. The time varying gyro bias, b, can be 
obtained by the following shaping filter : 


dz 

sin u sin i + — cos u 
at 

d i 

cos u sm i — — sin u 
dt 


(3.9) 


cosi + u. 


dh 

dt 


m, 


(3.12) 


where fj 2 is another white noise uncorrelated to fj\ . 

The noise characteristics of HRG under consideration for the ICESAT are given in 
Table 3.3. The components of Equation 3.11, such as rj \ , 60 , and b, can be derived from 
these values. 


3.3 Star Catalog 

3.3.1 Star Catalog for Star Field Simulation 

The star catalog is a fundamental part of the attitude determination process that uses 
measurement data obtained from any type of star sensor. For the star field simulation, 
a star catalog, originally developed to support a spacecraft equipped with a Ball CT-601 
star tracker, was obtained from the Smithsonian Astrophysical Observatory [56]. The star 
catalog contains 4853 stars and is a subset of Yale Bright Star Catalog. Two years of orbit 
simulation of ICESAT sampled the entire celestial sphere because of the complete rotation 
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Table 3.4: Probability of the number of stars in a FOV 



Whole Sky(4853) 

South Pole(391) 

North Pole(323) 

Stars 

Poisson 

Computer 

Poisson 

Computer 

Poisson 

Computer 


dist.(%) 

Simul.(%) 

dist.(%) 

Simul.(%) 

dist.(%) 

Simul.(%) 

< 2 

2.0 

3.9 

0.6 

2.8 

2.1 

1.1 

3 

3.8 

5.5 

1.4 

3.2 

3.9 

2.8 

4 

7.2 

8.6 

3.3 

6.0 

7.4 

8.4 

5 

10.8 

11.6 

5.9 

9.5 

11.0 

12.5 

> 6 

76.2 

70.4 

88.8 

78.5 

75.6 

75.2 


of node through 360°. The probabilities of a certain number of stars being in a FOV using 
the Poisson distribution and the computer simulation are presented in Table 3.4. From 
the main purpose of the GLAS mission, which measures ice sheet elevation over Greenland 
and Antarctica, the distribution of stars in both polar regions (above 60° and below —60° 
of declination) is more important than any other regions. Therefore the independent 
computations for polar regions were performed and presented in the table. The number 
of stars in each polar region is given in the parenthesis. If the distribution of stars were 
ideally uniform in the star catalog, the star tracker would observe seven or eight catalog 
stars in each image. However, the non-uniform distribution of stars in the star catalog 
(and on the celestial sphere) gives many different numbers of stars in each image as shown 
in Table 3.4. If there are more than six stars in the HD-1003 star tracker FOV, the built-in 
processor may select only six of them in terms of their brightness or relative positions. If 
there are one or two stars, the stars still can be identified using the direct match technique 
(DMT) which will be introduced in Chapter 4. However, the accuracy of the determined 
attitude is lower than that obtained from more stars. An approximate relation between 
the number of stars and the attitude determined by a deterministic method is [25] : 


O star 

P VNfov 


(3.13) 


where a p d is the accuracy of the pointing direction perpendicular to the detector plane, 
a s t ar is the uncertainty of the star position and Npov is the number of stars in the FOV. 

In the simulation, no cases occurred when no star was observed in the FOV. However, 
eclipses by the Sun or the Moon can produce periods with no star observations. The 
approximate ranges where the star tracker is adversely affected are approximately a 35° 
radius from the Sun and a 25° radius from the Moon [11]. For example, the 35° radius 
from the Sun may cause the maximum eclipse period of about 19 minutes. During the 
eclipses by the Sun or the Moon, the attitude determination based on the star tracker 
measurement is not available and then the attitude changes need to be detected by a set 
of gyros (e.g., IRU) until new star measurements are obtained in the FOV. However, if the 
time duration in which only gyro measurements are available is too long, the intrinsic bias 
drift of gyros will cause a significant deviation of the determined attitude from the truth. 
To reduce the influence of the Sun and the Moon, the ICESAT/GLAS PAD could use Ball 
CT-602 star trackers whose BDs are tilted in terms of the HD-1003. In the low declination 
the eclipsing due to the Sun and the Moon would be forecasted and pre-considered, but 
the low latitude region requires relatively relaxed attitude accuracy. For the ice-sheet 
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Table 3.5: Specifications of the Hipparcos Star Catalog [1] 


Measurement Satellite 

Hipparcos Satellite(ESA) 

Measurement Period 

1989.85-1993.21 

Number of entries 

118218 

Catalogue epoch 

J2000 

Reference system 

ICRS 

Date Published 

June 1997 

Astrometry (. Hp < 9 mag) 


Median precision of positions, J1991.25 

0.77/0.64 mas (RA/dec) 

Median precision of parallaxes 

0.97 mas 

Median precision of proper motions 

0.88/0.74 mas/yr (RA/dec) 

Estimated systematic errors 

<0.1 mas 

Photometry (Hp < 9m ag) 


Median photometric precision 

0.0015 mag 

Mean number of photometric observations 


per star 

110 

Number of entries variable 


or possibly variable 

11597 

Number of solved or suspected 


double/multiple systems 

23882 


measurements in polar regions, solar and lunar eclipsing will not raise a serious problem 
because the Sun and the Moon are located in low declination (—30° < 5 < 30°). Planets 
are also positioned on or near the ecliptic plane and their movements are well predicted. 
Thus, planetary obscuration will be dealt with in a similar manner to eclipses by the Sun 
and the Moon. 

3.3.2 Star Catalog for Real Application 

For the PAD, the recently completed Hipparcos Star Catalog [1], containing the most 
accurate astrometric and photometric star data compared with other star catalogs, will 
be used.* The median precision of position and brightness of the Hipparcos Star Catalog 
are 0.77 milliarcsecond and 0.0015 magnitude respectively. The important features of the 
Hipparcos Star Catalog are given in Table 3.5. 

* SKY2000 Master Catalog (Version 2) was used in practice. The Hipparcos and Tycho catalogs provide 
high-quality astrometric and photometric data for virtually all SKY2000 stars in Version 2 and later 
versions. The Tycho Catalog is also the output from ESA’s Hipparcos mission. 
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3.3.3 Corrections of Star Measurement 

Aberration 


Aberration is the apparent shift in the position of a star caused by the motion of the 
spacecraft. For Earth-orbiting spacecraft, the motion of the Earth around the Sun causes 
a maximum aberration of about 20 arcseconds. The motion of the spacecraft about the 
Earth at a 600 km altitude with a circular orbit accounts for about 6 arcseconds additional 
aberration. The aberration angle, j3, is computed from the spacecraft velocity relative to 
the Sun, v, by [65], 


/3 = — sin ( 

c 


(3.14) 


where c is the speed of light and 9 is the angular separation between v and the true star 
direction. Because we need information for all directions, the vector form of aberration 
equation is necessary. By using the nutation angles, the aberration vector, k, and true 
obliquity of the ecliptic, the unit vector direction to the star corrected for annual aberration , 
Sa, is [54] 


Sa = (1 + St ■ k)St - n, 


(3.15) 


where St is the unit vector to the star rotated into true-of-date coordinates from the unit 
star vector in mean equatorial coordinates of date. The aberration due to the spacecraft 
motion around the Earth is approximately computed from 


S = (1 — Sa ■ v/c)Sa + v/c, (3.16) 

where S is the unit vector to the apparent place of the star observed by the star sensor in 
the spacecraft. The aberration correction is applied right after the star identification. 


Proper Motion 

Many stars show continuous changes in position indicating a certain angular velocity. Such 
angular velocities are known as proper motions. The proper motion of an individual star 
may be as large as several arcseconds per year. Special considerations must be given to 
the stars which will show significant changes during the mission period due to large proper 
motions. 


Parallax 

The star catalog is usually created in the heliocentric inertial coordinate system. Since the 
Earth is moving around the Sun once a year, the direction of a star as seen from the Earth 
(and the spacecraft) is changing periodically and half of the changed angle is called parallax. 
For the GLAS mission, the corrections to the parallax are not required since the maximum 
parallaxes for the very few closest stars to the solar system are only 0.8 arcsecond and 
parallaxes of most stars are negligible. Furthermore, the attitude determination is based 
on all the stars in the FOV so that the parallax error on one or two stars will not affect 
the result significantly [25]. 
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Chapter 4 


STAR IDENTIFICATION 


An essential component of CCD star tracker data processing is the star identification. 
Star measurement data must be identified using the star information in the mission star 
catalog to determine exactly which stars the sensor is tracking. Star identification al- 
gorithms require appropriate parameter adjustment depending on sensor characteristics, 
noise environments, and the given star catalog. There are several existing star identifica- 
tion techniques [65]. In the first section of this chapter, we will discuss a pattern matching 
algorithm (PMA), which matches the angular distances between pairs of observed stars 
with those of cataloged stars. Since the CCD star tracker enables us to detect multi- 
ple stars simultaneously, it seems appropriate to choose the PMA as a star identification 
method. An advantage of PMA is that this method can be used when no a priori attitude 
(or star) information is obtainable or the quality of the a priori information is in doubt. 
The PMA developed for this research requires at least three stars at a measurement time, 
but the simulation sometimes showed that only one or two stars were observed in a star 
tracker FOV. The second section of this chapter will discuss the DMT, which identifies 
every measured star separately in the star catalog. Since the ICESAT /GLAS will stay in 
the simple nadir pointing attitude mode and will estimate the attitude with one arcsecond 
accuracy at the measurement time, the predicted attitude would be close to the true atti- 
tude. In other words, we have good prediction for star positions at the new measurement 
time and then the area to be searched for in order to find the matched star should be 
small. In this way, it is possible to identify most of the measured stars, even when one 
or two stars are observed by the star tracker. The DMT could be used as the auxiliary 
method to help the PMA, but it could also be used as stand-alone method if the BD would 
be known all the time with sufficiently good accuracy either from the attitude prediction 
or from real time on-board attitude determination. 

4.1 Pattern Matching Algorithm 

4.1.1 Sectioning of Star Catalog 

The star catalog for attitude determination contains at least several thousand stars. There- 
fore, a search for the matching stars in the star catalog that occur within an 8° x 8° FOV 
might take an excessive length of time. This also increases the possibility of misidentifica- 
tion since there would be a significant probability of encountering similar distributions of 
observed stars in the star catalog. A search of the entire catalog will cause more serious 
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problems in real time attitude determination where fast processing is required. However, 
if we have a best guess of the star tracker BD available from a previous estimation or 
from other (coarse) sensors, only a small area in the star catalog around the BD would 
need to be searched to identify the measured stars. This will greatly reduce the required 
search time and the probability of misidentification. To support this technique, the celes- 
tial sphere must be divided into many cells (or segments) in an orderly pattern so that we 
can find the matched stars in several cells surrounding the estimated BD. Catalog stars 
will be preassigned in those cells by their positions in the star catalog. 

As a first step to set up the divided cells in the celestial sphere, the locations of the 
cell centers are given by [59] 

S n = “ — cos -1 (£ n ) n = 0, 1,2.. ..JV 

a n j = — j _ o i 2 2n (4.1) 

3 2n + l J y ’ 

where 5 is declination, a is right ascension, and N determines the total number of cells 
and the size of each cell. The £ n is defined by 

77 7 r 

e„ = (-1)" cos(^ TI ) n = 0, 1, 2....N. (4.2) 

These equations divide the celestial sphere into N + 1 declination zones and (2n+l) equally 
spaced regions in each zone, yielding (N + l) 2 cells without overlapping. 

If the cell size becomes larger, more stars will be stored in one cell and a smaller 
number of cells will be made. In contrast, the small cell size will require access to more 
cells in order to find the matching stars. The choice of cell size (i.e., N) affects the star 
identification efficiency both in time and in storage. The area of the star tracker FOV 
must be considered to determine the optimal cell size. If N is 22, each cell covers an area 
slightly larger than an 8° x 8° celestial area which is the size of the FOV for typical star 
trackers (e.g., CT-601 and HD-1003). Assuming that the estimated BD of the star tracker 
is close to the true BD, Figure 4.1 shows that we only need to look at the area of 24° x 24°. 
If the star identification algorithm begins to search from the nearest stars to the estimated 
BD, the full search of 24° x 24° would not be necessary in most cases. 

The total number of cells is 529 for N = 22. Except for both polar regions, the cells 
are all trapezoidal shapes. After dividing the entire sky into 529 cells by the rule given 
in Equation 4.1 and 4.2, the 4853 stars in the star catalog were assigned into those cells 
by the star coordinates. Figure 4.2 and 4.3 show the cells near the north pole region and 
the south pole region, respectively. The number of stars in each cell is indicated within 
parentheses. The numbering of the cells starts from the north pole and it continues to the 
south pole. This alternative numbering between north and south hemispheres continues 
until it arrives at the celestial equator zone which contains 45 equally spaced cells with all 
the centers located at 2° declination. The maximum number of stars in a cell turns out 
to be 37, which occurs once. Three cells contain one star in a cell, however, there is no 
cell that has no stars. As we mentioned earlier in Section 3.3, either cells with excessive 
number of stars or with too few stars make the star identification troublesome and may 
degrade the accuracy of the attitude determination. It might be possible to develop more 
complicated codes to deal with dense and sparse regions with different criteria. 

To access the desired cells with the knowledge of the estimated BD (in a and 5), the 
location of the cell is determined by a pointer, n 2 + j. where n and j are obtained from 
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a FOV 



+ is the estimated boresight direction 


Figure 4.1: A cell pointed by an estimated boresight direction and nearby cells 


[59] 

n = 2[(| - <$)/A<5 + 0.5] (5 > 0°) 

= 2A + l-2[(^-#)/A5 + 0.5] (5<0°) (4.3) 

j = [a/Aa + 0.5], (4-4) 

where Act and A 5 are the width of each cell in a and 5. The symbol [ ] indicates the 
smallest integer greater than (or equal to) the number inside. 

To access the surrounding cells, we need to tabulate the numbers of surrounding cells 
for each cell number. Figure 4.1 shows the surrounding cells around a cell pointed by the 
BD. Once the cell number is given, whether it is a pointed cell by BD or a surrounding 
cell, the stars residing in the cell will be referred to one by one. 

4.1.2 Algorithm 

The position angles (<f> and A in Equation 3.6 and 3.7) of the observed stars in a FOV of a 
CCD star tracker are converted to the unit vectors by Equation 3.5. Two base stars can 
be selected arbitrarily. They might be the two brightest stars in the FOV to reduce the 
possibility of misidentification. The cosine of the angle, D ni , between a pair of measured 
stars is computed by 

D’; 2 = 5i • 5 2 , (4.5) 

where Si and S 2 are the unit vectors of the base stars expressed in the SCF. From a priori 
attitude knowledge, the cell pointed by the BD and the surrounding cells are known at 
the measurement time. The catalog stars in those cells are considered as candidate stars 
for identification. With the assumption that the estimate of the BD is close to the true 
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24 ( 12 ) 



Figure 4.2: Catalog cell pattern near the north pole. The number of stars in each cell is 
given inside parenthesis. 


BD, the candidate stars are paired by the distance from the estimated BD. The cosine of 
the angle between two paired catalog stars, i and j, is 

D? = S i -S j , (4.6) 

where S t and Sj are the unit vectors, defined in Equation 3.2, expressed in the CRF. Thus, 
the match is considered as found if the condition 

\D^-D^\<e (4.7) 

is satisfied, where e is an error window depending on both the star tracker measurement 
error and the star catalog position uncertainty. If e is too large, the possibility of a 
misidentification increases. If it is too small, a misidentification is unlikely, but the number 
of identified stars are apt to decrease. 

It is reasonable to say that the selected catalog stars, S) and Sj, likely match the 
measured pair of stars, Si and S 2 , when Equation 4.7 is satisfied. To resolve an inevitable 
180° ambiguity stemming from the angular distance comparison, the magnitude test is 
performed. This is processed in two steps. First, the two catalog stars i and j must have 
a magnitude difference greater than the magnitude error bound, /i. The /i is determined 
by combining the star tracker magnitude error and the catalog magnitude uncertainty. 
Second, the magnitude of the catalog star i must be close to that of the measured star 1 
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Figure 4.3: Catalog cell pattern near the south pole. The number of stars in each cell is 
given inside parenthesis. 


(or 2) within /i, while being different from the magnitude of catalog star 2 (or 1) more than 
/x. If these conditions are satisfied, the base stars, 1 and 2, are temporarily considered to 
be identified. 

Failure to resolve 180° ambiguity by the magnitude test will occur in two situations. 
First, the magnitudes of catalog stars may be unmatched with those of the measured stars 
in any combination. In this case, the catalog star pair, which satisfied Equation 4.7, are no 
longer considered as possible matches for the measured star pair. The catalog star pair are 
discarded and other pairs of candidate stars are examined sequentially. Second, when the 
brightness of the two catalog stars (or two measured stars) are too close, their separation 
with the parameter /j may not be obtainable. In this situation, the third identified star 
is required to solve the 180° ambiguity problem of the base stars. Even though the base 
stars pass the magnitude test, the third identified star is still necessary because there is 
relatively high probability of finding an invalid star pair match produced by too many 
catalog star combinations, the unregistered (on star catalog) background stars in the true 
sky, nearby space debris and/or ghost images [67] of the CCD star tracker. For this reason, 
the star identification is not initiated when only two stars (or less) are in a frame of the 
star tracker, and that frame is simply discarded. 

The third star is searched for in the remaining candidate stars, whether the 180° 
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ambiguity of base stars are cleared or not, until a star satisfying the following conditions 
is found : 


\D^-D^\ < e 

\Dt’ k -D%\ < £ , (4.8) 

where k and 3 are the indices for the third star in the catalog and in the measurement, 
respectively. If the magnitude test for base stars is successful with the matches of (i, 1) 
and (j, 2), Equation 4.8 is the condition which must be satisfied for the candidate star to 
be the third identified star. In contrast, if there remains a 180° ambiguity problem in the 
base stars, another chance exists to match the third star by exchanging the order of the 
base stars using the conditions : 


\Di k -D 2 J\ < e 

\D^ k -D^\ < e. (4.9) 

If the third star satisfies Equation 4.8 or Equation 4.9 and the measured magnitude of the 
third star is within the error bound of the corresponding catalog magnitude, all three stars 
are presumed to be identified. After the third star has been identified, the fourth star is 
matched in a similar way that was used for the third star. This process will be continued 
until all the measured stars (up to five or six stars) are checked with the candidate stars. If 
an observed star is not matched with any catalog star in this procedure, it is considered to 
be a false measurement and it is ignored. When only three stars are identified in a FOV of 
the CCD star tracker, there remains a chance of misidentification which will be discussed 
in the next section. If four or more stars are identified, the probability of misidentification 
is extremely small (negligible) primarily because the stars are searched in a limited number 
of cells with an orderly process and the estimated BD will be very close to the true BD. 

The values for e and n are critical for the star identification. Those values can be 
initially determined by the noise characteristics of the star tracker measurements and 
the star catalog uncertainty, and must be adjusted in the real data process after launch. 
The outline of the star identification algorithm developed in this section is illustrated in 
Figure 4.4. 

4.1.3 Misidentification 

Using the PMA, the measured stars can be absolutely identified or unidentified. For 
the latter case, the angular velocity detected by the gyros will enable prediction of the 
spacecraft attitude until identified star data are available. When only three stars are 
identified, there is a possibility of misidentification if they satisfy 

D 2 J = D 1 ^. (4.10) 

This is the condition that three identified stars form an isosceles triangle with the third 
star on or near the vertex in the symmetric line, as illustrated in Figure 4.5. When two 
base stars are located at the positions 1 and 2 and the third star is located at position 
3, there is little chance of a misidentification. However, the misidentification likely occurs 
when the third star is located on or near the symmetry line between the base stars. Then 
the three stars (1,2 and 3) form an isosceles triangle. If the distinctive magnitudes of 
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Figure 4.4: Outline of a star identification algorithm 
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Figure 4.5: Example of a star distribution in a CCD star tracker frame. An isosceles 
triangle is drawn by base stars (1 and 2) and a star at 3' (or 4'). 


the base stars are known by the magnitude test, the third star located at or near the 
symmetry line provides more confidence on the correct identification of the base stars. 
However, as long as the 180° ambiguity of the base stars is not resolved and if the third 
star is located at or near the symmetry line, insufficient information exists to distinguish 
the base stars. That measurement frame is simply discarded and rendered as unidentified 
at the expense of the loss of the possible correct matches. Since it is not likely to have two 
(for example, 3' and 4' in Figure 4.5) or more stars located on or near the symmetry line, 
the case of only three identified stars would need to be checked for the possible source of 
the misidentification. Therefore, procedures to detect the isosceles triangle with a 180° 
ambiguity must be added to the star identification algorithm. 

In addition to the isosceles triangle shape, there could be other special geometric 
configurations made by observed stars. For example, if there is a background star in 
position 4" (symmetric about the base line to the third star in position 3), it may be 
misidentified as a star 3 regardless that the star ( 4 ") is a registered star in the star catalog 
or a simple background star. We can imagine other special geometries formed by three, 
four or five stars which require more sophisticated star identification algorithms to avoid 
misidentification. 

There are other sources of misidentification that originate from non-geometrical rea- 
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sons. One example is space debris (or other satellites) near a star tracker BD. Since the 
star tracker does not discern these sources from a real star, debris may cause nrisidenti- 
fication. (Debris actually caused a more serious problem in the star tracker performance 
of the X-ray Timing Explorer (XTE) spacecraft [15]). Ghost images of the CT-601 have 
been reported as an identification error source [67]. The star identification algorithm must 
be sufficiently smart to consider ghost images as background stars that are unregistered 
in the mission star catalog. 

It is possible in the attitude estimation procedure to detect the measurement frame 
which has misidentified stars. The abnormal discontinuity of quaternions (or the com- 
puted angular velocity) obtained by Single Frame Attitude Determination (SFAD) meth- 
ods (Chapter 5) could be interpreted as an outbreak of the misidentification. Even though 
we suppose that the PMA could include the detection of various misidentification sources, 
the complete escape from the misidentification would not be easy. 

4.1.4 Simulation Results 

The star tracker measurements were simulated for four circular orbits with 94° inclina- 
tion using the algorithm developed in Section 3.1. The respective initial positions are 
0°, 45°, 90° and 135° right ascensions at the equator. The time interval between succes- 
sive FOVs is 0.1 second. The measurement noises were computed based on the lu values 
given in Section 3.1.3 for generating realistic simulation data. 

The PMA was applied to the simulated star measurement data and the results are 
summarized in Table 4.1. In the simulation, the possible maximum number of FOVs 
for each orbit is 57901 corresponding to a 10Hz measurement rate for one GLAS orbital 
period. Since the minimum number of stars necessary for the identification process is three, 
the number of frames containing at least three stars are shown in the second column of 
Table 4.1. The II in the third column is a parameter representing the range in which the 
third star is considered in a position to form an isosceles triangle with base stars. It is 
computed by 

U=\D 1 J-D 2 J\ (4.11) 

which means the distance difference between two sides of triangles other than the base 
line. In other words, if the third star, 3, is located at a position to make the right side 
of Equation 4.11 less than II, the three identified stars are considered to form an isosceles 
triangle. These stars are regarded as possible misidentification cases and are abandoned. 
The identification processing is resumed in order to find another group of stars which do 
not construct the isosceles triangle defined by Equation 4.11. 

Table 4.1 shows that all the star tracker measurement frames are reported as identified 
(indicated by 100%) if the isosceles triangle formed by the three stars is simply regarded 
as a correct identification (II = 0). However, the misidentification occurred in some frames 
as indicated in the last column of the Table 4.1, except for the case with 45° initial right 
ascension, a. For example, 0° initial a contains 167 misidentified frames, which is 0.3% of 
the reported identified frames. Computed results show that when the angle n is increased, 
the occurrence of the misidentification is decreased and finally disappears. There are more 
reductions in the total number of identified frames (see the seventh column in Table 4.1) 
than the number of misidentified ones (the eighth column). This clearly indicates that 
the elimination of the misidentification source results in the removal of some correctly 
identified frames. Table 4.1 also shows that the number of frames with five identified 
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Table 4.1: Star identification simulation results using the pattern matching algorithm. 
One orbital period of simulated data were processed. 


RA 

FOVs 

n 

Number of Identified Frames 

Mis. ID 

(a) 

(>3) 

(radian) 

3 stars 4 stars 

5 stars 

total(%) 

frames (%) 



0.0 

4418 

6074 

44487 

54979(100) 

167(0.3) 

0° 

54981 

2.0 x 10“ 4 

4399 

6074 

44487 

54960(99.96) 

167(0.3) 



1.0 x 10 -a 

4296 

6074 

44487 

54857(99.77) 

167(0.3) 



5.0 x 10 -a 

3477 

6074 

44654 

54205(98.59) 

0(0) 



0.0 

2729 

4186 

46704 

53619(100) 

0(0) 

45° 

53619 

2.0 x 10 -4 

2717 

4186 

46704 

53607(99.98) 

0(0) 



1.0 x 10“ a 

2655 

4186 

46704 

53545(99.86) 

0(0) 



5.0 x 10 -3 

2577 

4186 

46704 

53467(99.72) 

0(0) 



0.0 

2450 

5008 

49533 

56991(100) 

72(0.1) 

90° 

56991 

2.0 x 10" 4 

2441 

5008 

49537 

56986(99.99) 

68(0.1) 



1.0 x 10 -a 

2401 

5008 

49552 

56961(99.95) 

53(0.1) 



5.0 x 10 -a 

2190 

5008 

49605 

56803(99.67) 

0(0) 



0.0 

1666 

3821 

50394 

55881(100) 

115(0.2) 

135° 

55881 

2.0 x 10 -4 

1662 

3821 

50394 

55877(99.99) 

115(0.2) 



1.0 x 10 -a 

1642 

3821 

50394 

55857(99.96) 

115(0.2) 



5.0 x 10 -3 

1299 

3821 

50468 

55588(99.48) 

41(0.1) 



1.0 x 10“* 

1258 

3821 

50509 

55588(99.48) 

0(0) 


stars (the sixth column) increases as II increases. This means that the misidentification 
in some frames hinders the determination of correct star matches. The lowest probability 
of successfully identifying stars without misidentification was 98.59% in the simulation. 
Even though this result is undoubtedly promising, further star identification tests (and 
necessary improvements) should be carried out with more background stars in which the 
magnitude range covers up to eighth or ninth. 

4.2 Direct Match Technique (DMT) 

4.2.1 Algorithm 

The requirement of more than two stars as well as the relatively high possibility of the 
misidentification when only three stars are present in the star tracker FOV were pointed 
out as shortcomings of the PMA. The DMT enables us to identify the stars when the 
FOV contains only one or two stars. Once the GLAS attitude determination is initialized 
and determines attitude in the required accuracy, the DMT can identify almost all stars 
observed by the FOV. The prediction of attitude can be achieved by applying the kinematic 
equation of attitude quaternions, Equation 2.11. Alternatively, the quaternions obtained 
by the real time on-board attitude determination can be used. The accuracy requirement 
for real time is 10-20 arcseconds(l(j). Figure 4.6 illustrates the relations of the coordinate 
frames in the DMT. 

The star tracker observes up to six stars at a measurement time. If either the predicted 
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or the real time attitude is known, the star positions will be transformed from the SCF 
to the CRF (the observed CRF in Figure 4.6). If either the predicted or the real time 
attitude is sufficiently close to the true attitude, the small circular area centered by the 
observed star will be searched in the star catalog (Figure 4.7). By defining e as the error 
window radius of the circular area, the catalog star satisfying the equation 

(1(0, S) < e (4.12) 

will be considered as the matched star, where d(0 , S ) is the angular distance between O, 
the observation unit vector in the observed CRF, and S, the catalog star unit vector in 
the true CRF. The value for e is determined by the accuracy of the predicted or real time 
attitude. 

It is easy to see that the procedure described here would provide unique star identifi- 
cation if the density of star population is considered. Since the HD-1003 star tracker will 
observe stars brighter than magnitude 6, the number of stars in the corresponding star 
catalog is generally five or six thousand. In this case, the stellar density of the catalog is 
roughly less than 1 star per square degree. By selecting several tens of arcseconds as the 
error window radius, e, in Equation 4.12, we are able to identify almost all stars uniquely. 
If multiple stars are found in the error window radius, additional work has to be done with 
the comparison of either star magnitudes or the error sizes, d(0,S), in order to achieve 
the unique identification. If no star is found within the radius e, the observed star might 
be the non-star object or the predicted (or real time) attitude might have large error. 
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Figure 4.7: Catalog stars close to measured stars 


4.2.2 Catalog Star Search 

In the PMA, the area to search for the matched star should be larger than the star 
tracker FOV, 8° x 8°. It was required to divide the celestial sphere into many cells in 
order to look at the restricted region in stead of the whole sky. (see Section 4.1.1). In 
the DMT, a star can be identified by looking at the small circular area whose radius is 
usually tens of arcseconds with the observed star at the center. The method sectioning 
the star catalog into many cells could be applied to the DMT, but the size of each cell 
must be much smaller than that for the PMA in order to efficiently search the area of the 
several tens of arcseconds e radius (Equation 4.12). This would require excessive work to 
create and organize more than several thousand cells and access to each cell in the search 
process would be problematic when we consider the limited computation time. Therefore, 
the DMT developed for the GLAS PAD applies a new approach that does not use the 
sectioning of the star catalog described in Section 4.1.1. As the first step of the new 
approach, the catalog stars whose right ascensions, a , are sufficiently close to the observed 
stars will be selected as candidate stars for identification. Since the star tracker is able 
to track the stars, the most probable matched star is the one that was identified at the 
previous measurement time. By looking at the catalog stars in both the increasing and 
decreasing direction of catalog numbers, from previously identified number, k, it is very 
likely to find catalog stars whose right ascensions, a, are in several tens of arcseconds from 
that of the observed star, unless the star is new to the star tracker FOV. Therefore, the 
search area can be restricted to : 


(Xq — 6 <C Oi <C (X o V c, 


(4.13) 
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(a) for existing star 



(b) for new star 



Aa, > 4° 
Aa 2 > 8° 


Figure 4.8: Search for matched star in star catalog 


where cto is the observed star a (see Figure 4.8(a)) and e is the error window radius already 
introduced in Equation 4.12. 

For a star newly moving into the FOY (see the bottom of Figure 4.7), the e boundary 
doesn’t help find the match and we need to extend the search area. We refer to one of 
the previously identified star catalog numbers and need to search ~ 8° right ascension 
range both in increasing and decreasing k directions (Figure 4.8(b)). To reduce the search 
time, we may further restrict the search area within ~ ±4° right ascension range from 
the predicted BD. Equation 4.12 is applied to the stars residing in the selected a range in 
order to choose the correct match for the observed star. 

The search method described here performed well in the simulation. However, due to 
the nature of the longitude and latitude system for the celestial sphere, a large number 
of stars (sometimes all the stars in the star catalog) had to be looked at when the star 
tracker observed the vicinity of the polar regions. This search method has been improved 
by dividing the celestial sphere according to latitude zones as shown in Figure 4.9. 

This zoning of the star catalog is a reduced version of sectioning described in Sec- 
tion 4.1.1. To avoid looking up to nearby zones, each zone is overlapped with the neigh- 
boring zones. This could enable us to search only a rectangular area represented by thick 
lines in Figure 4.9(b) instead of the area indicated in 4.9(a). 

4.2.3 Simulation Results 

In the simulation of the DMT, we achieved almost 100% successful star identifications 
without any misidentification. Table 4.2 shows the results obtained by the DMT when the 
same simulation data as for the PMA were used. 
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(a) 


(b) 


Figure 4.9: Searching area for star identification. Latitude zones more than 4 are expected. 


Table 4.2: Star identification simulation results using the direct match technique. One 
orbital period of simulated data were processed. 


RA 

(«) 

Total Number of 
Observed Stars 

Number of 
Identified Stars 

Number of 
Misidentified Stars 

0 ° 

265738 

265735 (99.999%) 

0 

45° 

265003 

265003 (100.00%) 

0 

90° 

276872 

276867 (99.998%) 

0 

135° 

275330 

275232 (99.964%) 

0 
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Chapter 5 


ATTITUDE DETERMINATION 


5.1 Single Frame Attitude Determination 

In this section, several algorithms will be introduced to determine the spacecraft attitude 
from the observed unit vectors. Since these algorithms determine attitude from measure- 
ment data obtained by a frame of a CCD star tracker, they will be called Single Frame 
Attitude Determination (SFAD) methods in this document. 

Because the loss function J(A) in Wahba Problem (Equation 2.20) can be scaled with- 
out affecting the optimum attitude matrix, A op t, it is possible to normalize the weights 
such as 

n 

J> = 1 (5.1) 

i = 1 

with no loss of generality. Then the loss function can be easily manipulated to 

n 

J(A) = 1 -^caWfAVi 
1=1 

= 1 -9(A). (5.2) 

The loss function J(A) will have a minimum value when the gain function g(A) has a 
maximum value. The g(A) is given by 

n 

g(A) = J2 a i WjAVi 

i= 1 
n 

= ^aMAViWf] 

1=1 

= tr[AB T ], (5.3) 

where tr means trace of a matrix and the B-matrix or Attitude Profile Matrix , B , is, 

n 

B = Y, MUkf. (5.4) 

i = 1 

The B-matrix is very important since it contains all the information from the measurement 
and reference vectors. The construction of this matrix is the beginning of every solution 
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of the Wahba Problem in this section. Finally, the original Wahba Problem proposed in 
Equation 2.20 turns into the minimization problem of 

J{A) = 1 - tr(AB T ) (5.5) 

or the maximization problem of 

J'(A) = tr(AB T ). (5.6) 

There exist many SFAD algorithms. The simple description of these methods are 
given in the following sections, while the precise derivations of these can be found in the 
references. 


QUEST 


The QUEST algorithm is an improved version of the q-method developed by Davenport 
[65]. He has shown that the minimization of the loss function J(A) can be transformed 
into an eigenvalue problem of a 4 X 4 matrix, K : 


where 


K 


S-al Z 
Z T a 


S = B t + B 

n 

Z = J^diWiXVi 

i = 1 

<7 = trB. 


(5.7) 


(5.8) 


It was shown that components of the eigenvector corresponding to the biggest eigenvalue 
in K matrix are the attitude quaternions. Shuster’s QUEST algorithm [49] is a variation of 
the q-method in the last step of q-method. It avoids the necessity for solving the eigenvalue 
problem and saves computation time by using a simple Newton- Rhapson method. Through 
matrix manipulations, Shuster finally derived a fourth order characteristic equation for A 
as 


A 4 — a A 2 — /3A + 7 = 0, 


(5.9) 


where a, f3 and 7 are defined in terms of S, Z and a in Equation 5.8. Since A max is known 
to be very close to unity when the normalized weight is used, Equation 5.9 is solved for A 
with the Newton-Raphson method using A = 1 as an initial value. This is the point where 
the QUEST algorithm is preferable to the original q-method from a practical point of 
view. In the simulation performed for this study, only one or two iterations were required 
to solve Equation 5.9. Once A max is known, the corresponding optimal quaternions can 
be quickly computed from 


1 

4 * + 11-hi 2 



(5.10) 


where s and X are also obtained by combining the parameters in Equation 5.8. 
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QUEST Covariance Matrix 


A rigorous derivation of the error covariance matrix is given by Markley [28] for the unit 
vector measurement. The general assumptions for the errors in the unit vector measure- 
ment are as follows : 

• Because the vectors Vi and W t are constrained to be unit vectors, the error in any 
one of them must, to the first order, lie in the plane perpendicular to that vector. 

• The error vector has an axially symmetric distribution about the respective unit 
vector. 

• Since the attitude estimation error must be independent of the true attitude, the 
statistical value for the size of the error in the attitude quaternions is determined 
only by the star tracker measurement and the star catalog errors. In other words, the 
amount of attitude error must be the same, regardless of the spacecraft orientation, 
if we use exactly the same algorithm and the same measurement and catalog errors 
for the attitude determination. Therefore, a special assumption can be made such 
as 

Wi = Vi, (5.11) 

which means that the inertial and observation vectors are identical. Most of the 
time, Equation 5.11 is not true. But this assumption does not affect the values of 
error covariance matrix (or 5q). 

Thus, the QUEST covariance matrix is derived as [49] 

1 n 

Pm= 1°ltlhx3 - Y. a iWiWl]~\ (5.12) 

1=1 

where ai is a weight of the measurement shown in the Wahba problem. The total mea- 
surement error in a frame of the CCD star tracker, <r to4 , is 

n 

2 tot )~ 1 = ( 5 - 13 ) 

i = 1 

where a* is the sum of the measurement and the catalog errors and is related to the weights 
of each measurement by 

°i = otot/\fai- (5-14) 

Applying the small angle approximation of quaternions in Equation 2.4, the covariance 
matrix for attitude error angles becomes 

Pee = 4 P m . (5.15) 

The square roots of the diagonal components of Pee correspond to yaw, roll and pitch (la) 
errors, respectively. 

Singular Value Decomposition SYD) 
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It is well known that any matrix, B, can be decomposed as follows : 

B = USV T , (5.16) 

where U and V are 3x3 orthogonal matrices and 

S = diag(si,s 2 ,s 3 ) (5-17) 

with 

si > S2 > S 3 > 0. (5.18) 

The Equation 5.17 means that S is a diagonal matrix with diagonal elements si,S 2 and 
S 3 . They are called singular values of B and Equation 5.16 is called Singular Value 
Decomposition (SVD) [58]. 

Assuming that the matrix B in Equation 5.16 is the same as the B-matrix defined in 
Equation 5.4, Markley derived the attitude matrix, A, as [28] 

A = U[diag(l,l,d)]V T , (5.19) 

where 

d = ( det U)(det V ) = ±1. (5.20) 


Polar Decomposition (PD) 

A very different approach to the Wahba Problem can be given as follows [ 8 ] : The Euclidean 
norm of the general real matrix Q is defined by 

IIQII 2 = ^ Qij = tr (QQ T )’ (5-21) 

where ||Q|| is the matrix norm of Q and the sum is over all the matrix elements. From 
the assumed orthogonality of A and the properties of the trace, the matrix norm of the 
difference of two matrices, A and B, is 

\\A - B || 2 = tr[(A - B)(A - B) T } = 3 - 2 tr{AB T ) + \\B\\ 2 . (5.22) 

Because the orthogonal matrix A that maximizes tr(AB T ) minimizes this norm, the 
Wahba problem is also equivalent to the problem of finding the proper orthogonal matrix 
A that is closest to B in the Euclidean norm. 

The Polar Decomposition says that every real square matrix can be factored into [58] 

Q = PR, (5.23) 

where the P matrix is orthogonal and the R matrix is symmetric positive semidefinite. If 
Q is invertible, then R is positive definite. This P is the closest orthogonal matrix to the 
matrix Q. Assuming error-free observations W, (i = 1,2 ,..n), Equation 5.4 becomes [ 8 ] 

n n n 

B = J2 aiWiV? = J2 aiAViV? = A^2 = AR, (5.24) 

i = 1 i= 1 i = 1 
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where R is a symmetric matrix. Therefore the attitude matrix A can be found by 

A = BR~ l . (5.25) 

A is not necessarily orthogonal because of the error-free assumption for attitude measure- 
ment. 

Iteration Method (ITER) and Improved Polar Decomposition (IPD) 

From the fact that the attitude matrix A is the closest orthogonal matrix to the B-matrix , 
the iterative orthogonalization algorithm [8] 

A 0 = B 

A+i = \{AT T + Ai ) (5.26) 

can be applied until Ai + \ converges to the solution of the Wahba problem. 

The PD method can be improved by applying the orthogonalization algorithm in Equa- 
tion 5.26 [8], then 

A = \{B~ t R + BFT 1 ). (5.27) 

This is called the Improved Polar Decomposition. 

Fast Optimal Attitude Method (FOAM) 

Combining the singular values of the B-matrix (see Equation 5.17), the optimal attitude 
matrix can be rewritten as [29] 

A = [(« + \\B\\ 2 )B + A adj ( B T ) - BB T B] / C, (5.28) 

where the scalar coefficients k, A and ( are determined using the singular values of the 
B-matrix. To derive A, the Wahba’s loss function is set to 

J(A) = A 0 - tr(AB T ), (5.29) 

where a* is not normalized, 

n 

A 0 = (5.30) 

2=1 

The benefit of the FOAM over the SVD method is that the coefficients k, A and £ can also 
be computed without using the SVD. By matrix algebra, k and ( could be represented as 
functions of A and B, where A is the root of the following equation : 

0 = p(A) = (A 2 - ||B|| 2 ) 2 -8 XdetB-4 ||adjB|| 2 . (5.31) 

The largest root of p( A) minimizes the loss function J(A). Since Ao is close to A, we have 
a very good choice of the initial A and need to iterate only a few times. 

Modified Fast Optimal Attitude Method (MFOAM) 
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All methods except QUEST (and the q-method) discussed up to now in this section produce 
rotation matrices rather than quaternions. The quaternions can be extracted from the 
attitude matrix using Equation D.5. 

The MFOAM [30] is the variation of FOAM in order to compute quaternions directly. 
Inserting Ao into Equation 5.28, the approximately orthogonal matrix, Aq, is obtained by 

A 0 = M/C(\ 0 ,B), (5.32) 

where 

M = ^(Ao + \\B\\ 2 )B + A 0 adj B T - BB T B. (5.33) 

Then, the attitude matrix can be orthogonalized while extracting normalized quaternions 
q from Aq. The extracting equations are a modification of Equation D.5. 

5.2 Extended Kalman Filter 

5.2.1 Basic Algorithm 

The estimation algorithm for the EKF is based on Lefferts et al. [23] and Fisher et al. [16]. 
Instead of the direct estimation of the four element quaternions, q, the angle error vector, 
59 : 

59 x 1 I" 5qi 

59 = 59 y =2 5q 2 (5.34) 

_ 59 z _ _ 5q 3 _ 

is estimated with the gyro bias error vector which is defined by : 

A b = btrue ~ b, (5.35) 

where b true is the true gyro biases and b is the estimated gyro biases. 5q\,5q 2 and 5q% 
are the vector components of small angle (or error) quaternions. The scalar component of 
error quaternions becomes 1, to the first order, when a small angle rotation is applied to 
the normal quaternions. In consequence, the dynamic model for the angle error vector is 
given as : 

—59(t) = [uj(t)]59(t) + Ab + ?7i (5.36) 

^A b{t) = fj 2 {t). (5.37) 

where [cc(t)] is an antisymmetric matrix which is the matrix representation of the vector 
product : 

[< a\b = —a x b. (5.38) 

Now, the conventional linear vector-matrix differential equation is specified by 

= F(t)5x + G(t)fJ(t), (5.39) 
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where 


m 

G(t) 

V 


[^fc-l|fc-l] ^3x3 
03x3 03x3 

hx3 03x3 
03x3 hx3 

m 


The propagation equation for the state covariance matrix P(t ) is : 


Pk\k-i — ®k-lPk-l\k-l®k-l + Qk- 1 ) 


(5.40) 

(5.41) 

(5.42) 


(5.43) 


where &k—i is the state transition matrix from tk - i to tk and Qk- 1 is the process noise 
covariance matrix at time tk- 1 - The subscript k\k — 1 denotes a prediction value before 
it is updated at time tk and k\k means the updated value with new measurement at time 
tk- The state transition matrix, <h(t), is obtained by 

^$(t,t k -i) = F(t)$(t,t k -i) (5.44) 

with the initial condition 

^{tk—ljtk—i) = 76x6, (5.45) 

where /6x6 is a 6 x 6 identity matrix. For our system, the state-error transition matrix is 


$fc-i = ${tk,tk- i) 

©fc-i ^fc-i 
03x3 hx3 J ’ 

where 

^0(Mfc- 1) = [u]@(t,t k _ l) 

^'S?(t,t k -i) = [w]T(t,t fc _i) + 7 3 x3 

subject to 

©(ifc-i)ifc-i) = hx3 

'I'(fyc-lTfc-l) = 03x3, 


(5.46) 

(5.47) 

(5.48) 

(5.49) 

(5.50) 


where O 3 X 3 is a 3 x 3 zero matrix. Numerical integration, such as a Runge-Kutta method, 
is employed to solve for the transition matrix, <&(t,tk- 1 ). The process noise covariance 
matrix, Qk- 1 , is used as a tuning parameter so that the optimal value may be found in a 
search to achieve the best attitude determination. 
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5.2.2 Observation Model 

Two separate measurement models can be used for the EKF. 

Unit Vector Filter (UVF) Observation Model 

The measured star position is converted into a unit vector which is assumed to have the 
same noise characteristics for all three components regardless of the angle that the vector 
makes with the BD [44] [47]. The corresponding measurement covariance matrix is defined 
as : 

Rk = o%h*3, (5.51) 

where a is the variance of the star measurement noise on each component of the unit 
vector at time t k . The measurement residual, z k . is defined by : 

4 = W k - W k \k-i- (5.52) 

W k is the unit vector measurement in the OBF and it is given by : 

W k = A k V k + AWVk (5.53) 

where V k is the unit vector of the same star in the CRF and W k \k - 1 is the prediction from 

W k \k-i = A k \k-\Vki (5.54) 

where A W k is a white Gaussian noise. By representing the rotation as a matrix exponential 
[48] : 

A k = e ^A k \ k _i 

~ (I-[^1)^ fc | fe -i. (5.55) 

Inserting Equations 5.53, 5.54 and 5.55 into Equation 5.52 : 

4 ~ [W k \ k _i}59k + A W k - (5.56) 

The sensitivity matrix, H k , can be derived as : 

H k = ([W fc | fc _i]i0 3 x 3 ). (5-57) 

QUEST Observation Model 

The QUEST is a deterministic algorithm to compute the attitude quaternions from si- 
multaneous vector measurements of celestial objects [49]. The vector measurement in the 
QUEST is given by : 

Wk : i = A k V k ,i + AIU fc> j, (5.58) 

where the subscripts, i and k , are the measured star index and the measurement time, 
respectively. The QUEST algorithm furnishes attitude quaternions, p k , as well as the error 
covariance matrix, R k , which is given by : 

n 

R k l = - (■ A kVk,i)(AkV k ,i ) T ], (5.59) 

i=l a k,i 
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where o\ • is the variance of star i’s measurement noise at time t k and n is the total number 
of the measured stars in the star tracker. The square roots of the diagonal elements of 
Rk are angle errors of the QUEST attitude. In the EKF, the QUEST quaternions, pk, 
are combined with the predicted estimated quaternions, qk\k-h using to build the O-C 
residuals , Zk, by 

z k = 2p k 0 q~l_ v (5.60) 

where <8> is a quaternion composition analogous to an algebraic addition of two successive 
rotations [65]. qk\k-i is obtained by using the attitude propagation : 


— q(t) = ~n(uj(t)) q(t ) 


fi(w) is : 


n(uj) = 


life — 1 

in the 

previous measurement 

0 

u z 

-UJy 

Ux " 


-Uz 

0 


UJy 


UJy 

^ X 

0 

UJ Z 



-UJy 

-UJ Z 

0 _ 



^ 0 . 01 ) 
l . and 

(5.62) 


The 


0-C residuals is divided into vector and scalar components : 


Zk = 


Zk 

1 ’ 


(5.63) 


where the last component of Zk is, to the first order, 1, because Zk is expected to be a 
small amount of rotation. The observation-state equation thus is 

Zk = S6 k + Vk 

= H k Sx k + Vk, (5.64) 


where Vk is a discrete white Gaussian noise vector originating from the QUEST measure- 
ment error. The sensitivity matrix, Hk, is 


Hk = [-13x3:03x3]- (5.65) 

5.2.3 Kalman Filter Update 

The EKF update equations are : 

Kk = Pk\k-i H k {HkPk\k-i H k + -Rfc) -1 (5.66) 

5x k \k = KkZk (5.67) 

Pk\k = {I ~ KkH k )Pk\k-i{I ~ KkHk) T + KkRkKk (5.68) 

for both observation models. Note that 5x k \k-i = O 3 X 3 in the EKF. For the UVF obser- 
vation model, the update procedure is applied only when the measurement time changes. 
Finally, the new estimates at measurement time tk are 

qk\k = ^ k \k®qk\k-\ (5-69) 

t>k\ k = frfc|fc-i + Ab k \ k- (5.70) 
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5.3 Batch Method 


The normal equations for batch estimation are : 

5x k = {HlR^H k + P^)-\HlR; l z k + P k l 5x k ) (5.71) 

P k = (HlR- l H k + P~ l )-\ (5.72) 

where 5x k and P k are the estimated state and covariance matrix at epoch time t k while 
5x k and P k are the a priori state and covariance matrix at the same epoch time. Without 
any previous estimates, P k is specified with arbitrarily large values. The terms containing 
the measurement covariance matrix, R k in Equation 5.71 and 5.72 are computed by 

n 

H%Rl l H k = (5.73) 

1=1 
n 

HlR- 1 ^ = (5.74) 

1=1 

where Hi is the same as H k in Equation 5.65 : 


Hi — [73x3.03x3] • 


(5.75) 


Ri is the measurement noise covariance matrix at time 1, and z) is the O-C residual vector 
defined in Equation 5.63 as z k . The transition matrix Q(ti,t k ) is obtained by : 

■ • • $(t k+1 ,t k ), (5.76) 


by computing transition matrices between all successive measurement times. This equation 
allows us to propagate attitude information from tj. to t k considering noisy gyro data at 
each measurement time. Therefore, the attitude angle errors and gyro bias errors are 
determined by Equation 5.71 : 


5x k 


so k 

A b k 


(5.77) 


at the epoch time, t k . Finally, quaternions and gyro biases at t k are obtained by : 


Qk,m = 256 k <g> q k 

,m— 1 (5.78) 

b k ,m — 1 T A6fc, (5.79) 


where m is the iteration number. These are actually the same form as Equation 5.69 and 
5.70 for the EKF, however, the states ( q and b ) are renewed over the previous values at the 
epoch time in the batch estimation. The iterative procedure continues until the following 
condition is satisfied 

I qk,m - q k ,m- i|<e, (5.80) 

where e is a pre-set value. Once the quaternions and biases are estimated at the epoch, the 
quaternions in every measurement time can be simulated by using Equation 5.61 through 
the time interval for batch. It is assumed that the estimated gyro biases are constants 
over a batch interval. At the next epoch, these gyro biases are used as initial values for 
another iteration. 
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5.4 Simulation Results 

The CCD star tracker, Ball CT-601 or Raytheon HD-1003, has a capability of observing 
five or six stars simultaneously within the tracker FOV. Only a maximum of five stars 
were processed in our simulation. It is expected that the determined attitude quality will 
be slightly improved by increased number of stars in a FOV, whether it is noticeable or 
not. 

Figure 5.1 shows the angle errors of the QUEST quaternions obtained from one orbital 
period of star tracker simulation data. Errors computed by comparison of the determined 
attitude with the true attitude show the noise-like behavior (the center part). There 
are two lines in the upper and lower parts corresponding to the computed lcr values. 
Figure 5.2 shows the angle errors obtained from MFOAM with the same simulation data. 
The attitude errors (and also the computation times) obtained by two algorithms do not 
indicate noticeable differences. The uncertainty in each of the roll and pitch axes is about 

3.5 arcseconds(lcr). The yaw error which is the angle error about an axis perpendicular 
to the tracker FOV is around 70 arcseconds(ler). The estimation error in yaw is always 
much bigger than those in roll and pitch because the yaw axis (the zenith direction) is 
assumed to be aligned with the star tracker BD. For the pointing knowledge of the laser 
beam toward the nadir, it is required to know accurate angle rotations about the roll and 
pitch axes, but not about the yaw axis. 

The results of attitude errors of all the SFAD methods discussed earlier are compared 
in Table 5.1. The five methods, QUEST, SVD, ITER, FOAM and MFOAM, do not show 
noticeable differences in the la values. The approximate CPU times on a Hewlett Packard 
(HP) model 735/125 Unix machine, with 125 MHz processor, are compared in Table 5.2. 
Two algorithms, QUEST and MFOAM, show relatively fast performance. While the CPU 
time is important in real time attitude determination and control, the differences shown 
in Table 5.2 might be unimportant in post-processing. As long as there are no distinctive 
differences in accuracy, however, it would be better to choose the algorithm with the fastest 
performance. 

Figure 5.3 shows the attitude angle errors from the EKF with the QUEST observation 
model. For the attitude propagation between star measurements, simulated gyro data 
were used. Noise-like true attitude errors and bounds (la envelopes) show significant 
improvement comparing to QUEST solutions (Figure 5.1). When less than three(< 2) 
stars are observed, gyro data are used to predict the attitude. During the absence of star 
tracker data, the la envelope could not be drawn in Figure 5.3. The large covariance 
values were obtained by the restart of the EKF with the new star tracker measurements, 
but quickly the bound was diminished to the steady-state values. 

Figure 5.4 shows the attitude angle errors obtained by the EKF with the UVF obser- 
vation model. The batch estimation was performed with a 300 second batch interval and 
the results are shown in Figure 5.5. The batch intervals are marked by vertical dotted 
lines every 300 seconds in the figure. The estimated quaternions at the epoch time were 
propagated for 300 seconds (until the following epoch time) and the propagated quater- 
nions were used to determine the angle errors shown in the figure. It can be seen that roll 
and pitch errors from the batch method are similar to those of the EKF methods. 

Finally, Figure 5.6 shows the EKF with the QUEST observation model, but the EKF 
did not use gyro data at this time. Since angular rate information is required for quaternion 
propagation, the angular velocity vector, instead of gyro bias vector, was included in the 
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Figure 5.1: QUEST angle errors in arcsecond 
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Figure 5.2: MFOAM angle errors in arcsecond 
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Table 5.1: SFAD output statistics. Mean errors are the averages of the determined yaw, 
roll and pitch errors from the known true attitude over one orbital period. The las are 
the averages of the square root of the diagonal terms in the covariance matrix. 


Method 

Axis 

Mean Error 
(arcsec) 

la 

(arcsec) 


yaw 

0.23 

67.95 

QUEST 

roll 

-0.00 

3.40 


pitch 

-0.02 

3.45 


yaw 

-0.04 

67.95 

SVD 

roll 

0.01 

3.40 


pitch 

-0.01 

3.45 


yaw 

-0.39 

112.72 

PD 

roll 

-1.50 

5.52 


pitch 

-0.27 

5.71 


yaw 

0.55 

67.95 

IPD 

roll 

0.04 

3.40 


pitch 

0.03 

3.45 


yaw 

0.11 

77.15 

ITER 

roll 

0.00 

3.86 


pitch 

-0.01 

3.92 


yaw 

-0.00 

68.17 

FOAM 

roll 

0.00 

3.41 


pitch 

-0.00 

3.48 


yaw 

0.00 

68.17 

MFOAM 

roll 

0.00 

3.41 


pitch 

-0.01 

3.48 


Table 5.2: SFAD computation time obtained by time command in HP 735/125 Unix 
machine with 125 MHz processor for one orbital period (~ 5800 seconds). 


Method 

user(sec) 

system(sec) 

real(min) 

QUEST 

92.8 

23.7 

1.58 

SVD 

99.0 

23.5 

2.03 

PD1 

93.4 

23.6 

1.58 

PD2 

94.3 

23.9 

2.01 

ITER 

107.9 

23.2 

2.13 

FOAM 

94.1 

23.4 

2.01 

MFOAM 

92.7 

23.4 

1.57 
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Figure 5.3: Angle errors from the EKF with QUEST observation model in arcsecond 
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Figure 5.4: Angle errors from the EKF with UVF observation model in arcsecond 


60 










0 2000 4000 6000 

time(sec) 

Figure 5.5: Angle errors from batch estimation in arcsecond 
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Figure 5.6: Angle errors from the EKF without using gyro data in arcsecond 
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Table 5.3: Computation time comparison 


Method 

Obs. Model 


user (sec) 

system(sec) 

real(min) 

EKF, 

QUEST, 

gyro 

149.4 

26.8 

3.00 

EKF, 

QUEST, 

non-gyro 

143.9 

24.7 

2.51 

EKF, 

UVF, 

gyro 

341.5 

10.1 

6.00 

BATCH, 

QUEST, 

gyro 

451.2 

92.0 

9.08 


state vector. Thus the equation of motions are now : 

—89(t) = [u;(i)]<50(i) + AcJ (t) (5.81) 

^Aw(t) = (5.82) 

where fj(t) is white Gaussian noise on the angular acceleration. Most of the remaining 
equations are similar to the previous ones for the EKF using the QUEST observation model 
except for some modifications caused by the decrease in the number of noise parameters. 
The simulation results shown in Figure 5.3 and 5.6 might raise a question about the 
usefulness of the gyro data. However, it should be noted that the simulation did not 

include high frequency spacecraft jitter that may not be measured by the star camera, 

but could be measured by the gyros. Until the high frequency jitter is considered in 
the simulation (or in actual data process), it cannot be concluded that the EKFs, case 
using gyro data and case without using gyro data, produce the comparable results for 
attitude determination. The approximate computation times, required to process one 
orbital period of data on a HP model 735/125 Unix machine, with 125 MHz processor, 
are shown in Table 5.3. Both EKFs with the QUEST observation model perform the 
computation faster than the EKF with UVF observation model. The batch method, as 
expected, is slower than any EKF. 

Figure 5.7 shows the relation between true estimation errors and a la bound obtained 
from the state error covariance matrix. The numbers on the line are associated with the 
tuning parameter, the process noise covariance matrix. As the number increases (12, 13 
• • • ), the magnitude of the process noise covariance matrix decreases. The EKF diverges 
when the process noise covariance matrix is too small. 

Figure 5.8 shows the relation between true estimation errors and root mean squares 
( rms ) of one orbital period EKF estimation. The trend of rms coincides with that of true 
errors. This particular case indicates the possibility of using rms values in the filter tuning 
of the actual data where the true attitude is not known. 
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Chapter 6 

POINTING DETERMINATION 


Figure 6.1 shows a conceptual sketch of the Stellar Reference System (SRS) which is being 
developed by Goddard Space Flight Center to measure the pointing angle of the GLAS 
laser beam to the required accuracy (less than 1.5 arcsecond, la), relative to the star 
field. To reduce the alignment flexibility originating from multiple error sources, such as 
spacecraft jitter and temperature variation, all the components of the SRS will be at- 
tached to the GLAS optical bench. The SRS consists of a Laser Reference Camera (LRC) 
and the Attitude Determination System. The LRC contains optics, called Lateral Transfer 
Retroreflectors (LTR), which are designed to extract and redirect onto the Laser Reference 
Telescope (LRT) both a portion of the outgoing laser pulse and a beam from a Collimated 
Reference Source (CRS). The CRS is rigidly attached to a side of the star tracker and pro- 
vides the alignment information of the star tracker. The LRT also observes a star (LRC 
star) which is brighter than 7.5 visual magnitude. All of the images observed in the LRT 
are projected onto a 0.5° x 0.5° FOV Laser Reference Sensor (LRS). After locating the 
LRC star on the star field observed by the star tracker, other illuminated images in the 
LRS will subsequently be transformed to the star tracker star field in order to determine 
the laser altimeter pointing direction in the SCF or the OBF. When the LRC star is not 
available, the CRS images projected to the LRS will be monitored to check the stability of 
the whole system on the optical bench while the rigid body assumption is being applied. 
While the GLAS laser is fired at a 40 Hz rate, the LRS sees the laser at 10 Hz. The 
shot-to-shot direction of the laser beam is recorded by the Laser Profiling Array (LPA) 
with a frequency of 40 Hz. The centroids from 40 Hz laser images observed by the LPA 
must be located in the LRS frame and consequently in the SCF. Ultimately, the proper 
coordinate transformations would give us the 40 Hz laser beam pointing direction in both 
the TRF and the CRF. The details of the SRS instruments, such as mechanical require- 
ments, design concepts, construction materials and statistical quantities, are provided in 
the GLAS Critical Design Review [60]. 

6.1 POINTING DETERMINATION ALGORITHM 

Once a star observed in either the LRS or by the star tracker is identified, the star appears 
as a static point source in the CRF (see Star Catalog in Figure 6.2). The identified 
tracker stars will allow us to determine the attitude, A(t), at time t by applying attitude 
determination algorithms. The attitude matrix transforms the star coordinates in the 
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Figure 6.1: Stellar Reference System 
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LRS Star Identification 



Star Identification 
(Attitude Determination) 


Figure 6.2: Pointing determination using the Stellar Reference System 

CRF into those in the instrument OBF. If the alignment matrix, S(t), which transforms 
the star coordinates in the OBF into those in the SCF, is temporarily presumed constant, 
the combined matrix, M(t) : 

M(t ) = S(t)A(t) (6.1) 

describes the coordinate transformation from the CRF to the SCF. Once the matrix, M(t), 
is determined at every measurement time t, the LRC star identified in the CRF can be 
transformed to a dotted line at the star tracker FOV (see Star Tracker in Figure 6.2). The 
tracks of the LRC star both in the star tracker FOV and in the LRS FOV are input data 
to the least squares method which will be introduced in the following subsection. If the 
star in the LRC is brighter than magnitude 6, it will also be seen in the star tracker so 
that the projected LRC star positions on the star tracker FOV can be verified/calibrated 
with the observed star by the star tracker itself.* 

6.1.1 Estimation of the LRS frame in terms of the SCF 

Figure 6.3 describes the overlap of the LRS FOV on the star tracker FOV, but the LRS 
FOV is oversized to well visualize the relation between the two FOVs. Ideally, the two 
frames share the common origin. In the real situation, there might be an offset between the 
two origins. This offset is represented as the coordinate of the LRS origin, (x 0 // se t> Hoff set), 
in terms of the SCF. In addition, the LRS FOV might rotate by angle 8 against the SCF. 
Finally, the magnitude ratio of two FOVs, F, might change if focal lengths vary. Once we 
know these offset parameters, (x 0 ff set , y 0 ffset, 6, E), any point in the LRS FOV can be 
projected into the star tracker FOV using the following equation : 

X star lrs 
y star lrs 

where ( x s tar L Rsi Ustar L Rs ) the coordinate of the LRC star in the LRS frame and 
(x s tar ST , V star S T ) * s the coordinate of the same star in the SCF. Considering the velocity 

‘Actual LRS performance was degraded because stars were only measured on the night side of the 
orbit. Laser pointing determination has been done using an alternative algorithm instead of the method 
described in this chapter. The new algorithm is described in Appendix B. 


= F 

cos 8 sin# 


% star st 

% offset 

— sin 8 cos 8 


V star st 

— Voffset 


(6.2) 
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Figure 6.3: Estimation parameters for coordinate transformation between LRS frame and 
SCF 

of the GLAS/ICESAT and the size of the LRS FOV, we will have approximately 80 points 
(i.e. 8 seconds) of the consecutive LRC star measurements in both FOVs. With these 
data, the offset parameters can be estimated by applying the least squares method [58] to 
the Equation 6.2. Finally, Equation 6.2 is used to locate the laser point in the star tracker 
FOV by inserting the estimated parameters, x 0 // se t, y of fset, 6 an d F : 

x laserLRS 

yiaser LRS 

where ( xi aserLRS , yi aS er LRS ) are the coordinates of the LRS laser in the LRS frame and 
(, xiasersri yiaser ST ) are the coordinates of the same laser point located in the SCF. Equa- 
tion 6.3 is also used for the projection of the CRS point. After the laser points observed 
in the LRS FOV are located in the SCF, direction cosines (i.e. unit vector for the laser 
pointing direction) can be restored by considering the focal length of the star tracker. 
The construction of the laser pointing vector will be made once the laser travel time (i.e. 
range) is determined. The inverse of the matrix M(t), M(t) -1 , should be applied to 
transform the laser coordinates from the SCF to the CRF. As discussed at the end of 
Section 1.1, the laser spot position vector will be obtained by appropriately combining the 
laser pointing vector with the laser pulse travel time and the POD position vector. The 
spot location will be produced in the Earth-fixed axes, specifically ITRF by the proper 
coordinate transformation [42]. 

A simulation result for GLAS laser pointing determination is shown in Figure 6.4. The 
LRC star measurement data were created from a star catalog which contains 32511 stars 
ranging down to magnitude 7.5. The LRC star measurement noise was in the range of 0.2 
to 0.4 arcsecond depending on the star brightness [60]. In addition to the random noise in 
the overall SRS, the true shot-to-shot variation of the laser beam direction was given as 
0.3 arcsecond (lcr). The pointing errors in the figure are the absolute differences between 
the true and the estimated pointing directions. For the most part, the pointing errors were 
below 1.5 arcseconds, satisfying the accuracy requirement. 


= F 

cos 9 sin# 


%lasersT ^ offset 


— sin 0 cos 6 


yiasersT Voffset 
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Laser Pointing Error (10 Hz) 



time(sec) 

Figure 6.4: Pointing determination simulation result (10Hz) 


6.1.2 Absence of the LRC Star 

As shown in Figure 6.4, a half degree LRS FOV will occasionally observe stars. The time 
duration before a new star comes in to the LRS appears to be within a few minutes in 
the figure. To determine the maximum or average time between LRC star observations, 
simulations covering the whole sky might be necessary. The probability computation of the 
star tracker’s star observations, from both theory and simulation, indicated that the stars 
in our universe do not deviate considerably from the random distribution (see Table 3.4). 
Using the Poisson distribution [63], the computation shows that, during an orbital period 
of the GLAS, a star will be observed in the LRS more than 100 occasions, meaning at least 
a star per minute. Therefore, the LRS would almost always observe a new star within 
several minutes after the dropout of an old one, especially in the polar regions where no 
solar or lunar blocking would occur. Expecting that the line of sight stability for the 
GLAS optical bench would be less than 0.1 arcsecond(l(j) during one orbital period [34] 
[43], we conclude that the offset parameters obtained during a LRC star observation can 
be used after the star leaves the LRS FOV. When a new star comes into the LRS, the 
parameters will be updated and used until another new star is observed. If a noticable 
change of instrument alignment between star tracker and LRS occurs, it would be detected 
by monitoring the CRS points since the LRS also observes the CRS images projected from 
the star tracker CRS. Figure 6.5 shows the simulation results after filling up intervals 
between LRS observations using constant offset parameters before a new star. 

6.1.3 LPA measurement projection into LRS frame 

While the LRS sees star, laser and CRS points at a 10Hz rate, the LPA observes the laser 
beam at a 40Hz rate. The position of the laser in the LPA needs to be located in the LRS 
in order to determine the 40Hz laser beam pointing direction in the CRF. Although the 
LPA consists of an array of 80 by 80 pixels, only 20 by 20 array containing the laser point 
near its center will be read. An example of the laser profile is shown in Figure 6.6. 
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Figure 6.5: Pointing determination simulation result (10Hz, filled) 
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Figure 6.6: Laser profile observed at the LPA 
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A method to find the centroid (. Xi t LPA,Vi,LPA ) of the laser image is : 


x i,LPA 


Vi, LPA 


V^20 v^20 t 

/ 777/ Im,n 

V^20 v^20 t 

/ v i / > 1 ^ ^ Im,n 

Y^20 \~>20 t 

/ v ] ^ 2-/1 ^ -*m,n 

V^20 v^20 r 

/ , 1 / , 1 ^ Lfl lm,n 


(6.4) 


where I m ^ n is the laser intensity of m row and n column. 

Pre-launch calibration will provide the formula to project the LPA laser point to the 
LRS frame. However, the relative position of the LPA frame with respect to the LRS 
frame might change slowly due to launch load, temperature variation or instrument aging. 
The center of the LPA laser image needs to be located in the LRS FOV after launch. Since 
both the LRS and the LPA observe the same laser shots at a 10Hz rate (see Figure 6.1), 
these common measurements provide a way to project the LPA image center to the LRS 
coordinate frame. Measurements of laser in both LRS and LPA frames will contain both 
shot-to-shot noise and measurement noise. The mean position, ( xlpa , Ulpa ) is computed 
by: 


XLPA 

ULPA 


2^1 %i,LPA 

N 

sr~^N 

1^1 Vi, LPA 

N 


(6.5) 


where N is the number of total measurements, i is an index for the individual measurement 
and (x^lpa, Hi, lpa ) is described in Equation 6.4. Similarly, the mean position of LRS 
laser is calculated by : 


xlrs 

Vlrs 


l_j\ %i,LRS 
N 

2^1 Vi, LRS 
N 


(6.6) 


Finally, an individual LPA centroid, ( x i)L PA on LRS , Ui,LPA on LRs), is projected into the LPA 
frame by : 


Xi,LPA on LRS = K{xi,LPA - xlpa ) + xlrs 

Ui,LPA on LRS = K(y it LPA - Ulpa) + Vlrs, (6.7) 

where K is the ratio of size difference between the two frames. The projection procedure 
is illustrated in Figure 6.7 as a linear transformation of a coordinate system whose origin 
is located at the mean position of the LPA centroids. 

The laser point projection represented by Equation 6.7 is based on the assumption that 
the coordinate axes of the LPA are parallel to those of the LRS. The rotation angle of the 
LPA frame, 8, with respect to the LRS frame is illustrated in Figure 6.8 where the mean 

position of the LPA centroids is moved to the origin of the LRS frame. If we consider 

the rotation angle 9, the observed or uncorrected LPA laser point, ( Xi tUncar , y,; )Uncor ), will 
become the corrected one, (xj jCor , yi yCor ), by : 

%i,cor 
Vi, cor 



cos 9 

sin# 


%i,uncor 


— sin$ 

cos 9 


Vi,uncor 
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(6.9) 


Errors of LPA position, when the rotation is not considered, are computed by : 

A n = (Axf + Ay,j 2 )5, 

where 

A Xi — -^Lcor Xi,wncor 
A Vi — Vi, cor Vi.uncor • (6.10) 

Equation 6.8 implies that the size of A r* is determined by both 6 and 5ri, where 5ri is 
defined by : 

( x i,uncor 4"~ Vi, uncord 2 • (6.11) 

As we mentioned before, the uncertainty of 5ri on LPA measurements consists of 
two components: measurement noise and shot-to-shot variation. The peak shot-to-shot 
variation is expected to be one arcsecond [50]. The magnitude of the measurement noise 
would be less than one arcsecond (la) [32], The optical bench vibration and instrument 
distortion are expected to be 0.45 arcsecond in total (1 <t). The systematic component for 
those errors is not clearly quantified, but the long term line of sight instabilities of the 
instrument optical benches observed or required in other missions [20] [34] [43] lead us to 
believe that the value should remain within 0.5° (3 a) during 5 years in orbit even in the 
worst case. Inserting values 

5ri < 4 arcseconds(3a ) 

6 <C 0.5° (= 1800 arcseconds ) , 

into Equation 6.8, 6.9 and 6.10, we have 

A n = 0.06 arcsecond. 

This is clearly negligible when we consider other error sources in the SRS [32]. Therefore, 
Equation 6.7 would be sufficient to project the LPA centroids into LRS frame without 
worrying about the effect of the relatively rotated coordinate axes between two frames. 

6.1.4 40 Hz LPA measurements 

Equation 6.7 should apply to 40 Hz LPA measurement data even though the mean position 
of the LRS laser image will be obtained from 10 Hz measurement data. In order to use 
Equation 6.1 for transforming all of the laser points from the SCF to the CRF, the attitude, 
A(t), needs to be known at a 40 Hz rate. Since the attitude, A(t), is determined from 
star tracker and gyro data whose measurement rates are only 10 Hz, a linear interpolation 
is required to get the instrument attitude at 40 Hz. This interpolation can be justified 
only if the LPA frame is relatively fixed with respect to the LRS frame during the interval 
in which the interpolation is applied (i.e. 1 second). If the instrument or optical bench 
has significant high frequency jitter, the interpolated attitude might produce unacceptable 
pointing error, then it might be necessary to use gyros capable of 20 Hz or 40 Hz rate 
measurements. 

Figure 6.9 illustrates the overall measurement relations between star tracker, LRS and 
LPA. The final products of the PAD are laser pointing unit vectors in both ICRF and 
ITRF. 
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Figure 6.9: SRS overall procedure 

6.1.5 Collimated Reference Source (CRS) 

The CRS, which will be mounted kinematically to the reference surface on the star tracker, 
will maintain true pointing direction and collimation of its own laser pulse. The stability 
requirement is 0.2 arcseconds (lc) over a ±0.5°C gradient [60]. From the error sources 
in the relay of the CRS beam including LRT distortion and LRC centroid error, the 
measurement error in the LRC is expected to be similar to that of the LRC star. The 
monitoring of the CRS beam image over time will allow us to monitor the SRS structure 
stability and to verify the image quality in the LRC. Furthermore, during the loss of star 
tracker data due to, for example, solar and lunar eclipsing - the maximum duration is about 
20 minutes (see Section 3.3) - the offset parameters in Equation 6.2 may be presumed as 
constants if the CRS image remains in the stable position without any systematic change. 
The angular diameters of the Sun and the Moon are both approximately 0.5°, which is 
the size of the LRC FOV. As long as the LRC is not under the direct illumination of the 
Sun or the Moon, the LRC is expected to provide centroids of both CRS and laser images 
at the 10 Hz rate. By applying the algorithm depicted in the previous sections, the 40 Hz 
laser pointing determination will be available with the LPA laser images. 

6.2 Spacecraft Velocity Influence on Laser Pointing 

In Section 3.3.3, a simplified algorithm, to correct the stellar aberration, was shown. The 
PAD of the spacecraft is determined within the CRF and must be corrected for this effect, 
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stellar aberration being the result of the spacecraft (and the star tracker) motion with 
respect to the relatively fixed stars. In like manner, the motion of the spacecraft with 
respect to a point on the surface of the Earth will cause a similar effect that can be 
tentatively called pointing aberration. Figure 6.10 illustrates both stellar and pointing 
aberrations. Since the aberration is a function of relative velocity, and not the distance, 
between the observer and the source, the same algorithm can be applied. For pointing 
aberration, the applicable velocity is the relative velocity of the Earth with respect to the 
spacecraft. The magnitude of the pointing aberration can be calculated by Equation 3.14 

(3 P = — sin 0 
c 

I ^s/cl 

C 

= 5.2 arcseconds 

where the 9 is always close to 90° because of the near-nadir pointing direction of the laser. 

The pointing aberration will shift the angle of incidence from the laser direction that 
will be determined by the PAD processing. In other words, the laser spot will not be 
located at the path of the laser determined by the PAD. Instead, the laser will hit the 
ground at the path of the shifted direction. Without applying the effect of the pointing 
aberration the laser pointing direction determined from the GLAS instrument data will 
miss the actual laser spot by about 15 meters corresponding to 5.2 arcseconds. 
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Figure 6.10: Two aberrations due to the motion of the spacecraft 
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6.3 Systematic Errors 

6.3.1 Systematic Errors in the Stellar Reference System 

The pointing determination algorithm described in the previous section applied the inverse 
of the matrix M(t), M(i) _1 , to transform the estimated laser pointing direction from the 
SCF to the CRF. The laser pointing direction was determined in the SCF by referencing 
the LRC star trajectory transformed from the CRF to the SCF using the matrix M(t). 
Since we are ultimately interested in the laser beam direction in the CRF, the systematic 
errors embedded in the matrix M(t) might be canceled out by two way (i.e. direct and 
inverse) implication of the matrix M (t) to and from the SCF at a measurement time t. The 
matrix M(t ) is the combination of the attitude matrix, A(t), and the alignment matrix, 
S(t). The conceivable systematic error sources for the attitude matrix A(t ) are the star 
tracker CCD array pixel biases, gyro non-orthogonality and aberration correction errors. 
The alignment matrix S(t) can contain the systematic error due to the excursion of the 
star tracker BD. 

The excursion of the BD was added to the simulated star tracker data. The laser 
pointing direction was determined by the process introduced in earlier section. Figure 6.11 
shows the computed laser pointing errors (i.e. true laser pointing direction minus estimated 
laser pointing direction) for 600 seconds with various amount of the BD error. The first 
plot (a) is the result when no BD error was included in star tracker data. The average of 
the pointing error during the LRC star presence was 0.49 arcsecond. The addition of the 
200 arcseconds error on the star tracker BD did not change the pointing determination 
error. The average of the pointing error was again 0.49 arcsecond in the second plot 
(b). When the BD error becomes 1 degree, the average of the pointing error increased to 
0.84 arcsecond as seen in the third plot (c). With a BD error of 5 degrees, the pointing 
estimation no longer produces acceptable results as shown in the fourth plot (d). A similar 
simulation can be applied for the effect of the systematic errors included in the attitude 
matrix A(t). However, different approaches for the attitude systematic errors are discussed 
in the following section. 

6.3.2 Systematic Errors in Attitude Determination 

A Ball CT-601 star tracker was used for the high accuracy attitude determination of 
the Midcourse Space Experiment (MSX) mission [19]. The highest accuracy that the 
MSX had to achieve for the instrument pointing stability was 2 arcseconds (1 <t). While 
the pre-launch simulation showed the fulfillment of the attitude determination/pointing 
requirement, the MSX could not determine the attitude to the required accuracy from the 
actual star tracker and gyro data. In the attitude determination simulation before launch, 
some fundamental assumptions were made. Eventually, the assumptions turned out to 
be inadequate for the actual data processing. One of these assumptions concerned the 
characteristics of the star tracker measurement errors. While these errors were presumed 
to be adequately modeled as white noise in the simulation, the individual pixel of the 
star tracker CCD array included its own bias, eventually, resulting in the error of the star 
image centroid. Since the effect of the systematic error can not be removed by the usual 
statistical estimation procedures, the attitude determined by the white noise assumption 
could not satisfy the accuracy requirement. 
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Figure 6.11: Effect of the BD error in the Stellar Reference System 


While the MSX is completely versatile, with science scenarios varying from inertial 
stares to whole sky scans to tracking satellites, ballistic missiles and their surrogates, 
and covering about anything in between, the GLAS/ICESAT will maintain one attitude 
scenario : altimeter laser pointing toward the near-nadir direction in near-circular and 
near-polar orbit. A star would come in and move across the 8° x 8° FOV with a relatively 
constant rate for about 4 minutes until it exits. By considering 512 by 512 pixels in the 
CCD array, a star will stay on a pixel about 0.25 seconds. Unlike the case of the inertial 
staring attitude where one pixel sees a star for long duration, the pixel’s local bias will 
affect the star measurement in two or three snaps, and then the star moves into the next 
pixel that would be characterized by its own size and direction for bias. Moreover, the 
attitude determination algorithms will utilize multiple stars which will be located at all 
different pixels in the tracker CCD array at a measurement time. In this situation, the 
change of pixel biases of five or six star observations should be considered as random and 
the residual effect of the systematic error sources would be negligible. 

By using the bias correction tables within the post facto attitude determination soft- 
ware, the Roentgen Satellite (ROSAT) mission corrected the star tracker cameras’ mea- 
surements in pointing mode [37]. For the GLAS, it is expected that the manufacturer 
will provide the built-in bias calibration function of the star tracker [50] to remove both 
the individual pixel biases and the large distortion of the CCD array before the data are 
actually taken. This self-calibration ability of the star tracker will, at least, reduce the 
size of the bias errors significantly, if not remove them completely. 

The SRS ability to cancel out the measurement systematic error was also discussed 
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in the previous section. The laser pointing determination algorithm using the specially 
designed SRS will avoid the systematic errors embedded in the alignment matrix S(t), as 
well as in the attitude matrix A(t). When we consider only instrument pointing deter- 
mination, the SRS should be able to remove the effects of many systematic error sources 
even at the inertially staring attitude mode of the spacecraft. 

Besides the white noise model for the star tracker measurements, another inadequate 
assumption for the MSX attitude determination system was the obtainability of the un- 
interrupted star tracker attitude data at its nominal update rate of w 9 Hz. Extensive 
analytic studies indicated fairly strongly that glint from the Sun is a significant problem, 
possibly because of unanticipated reflections off the SPIRIT III (one of the MSX’s instru- 
ment) shield [19]. The unexpected glint would require the tracker integration time to be 
longer than anticipated, which will result in the large image smearing and interrupted data 
rate. The Kalman filters and the batch algorithm introduced in Chapter 2 are capable 
of estimating the gyro biases as well as the spacecraft attitude quaternions. The high 
frequency attitude update due to 10 Hz star tracker measurement rate will make the gyro 
biases include not only true bias but also other effects such as gyro non-orthogonality, 
gravity or magnetic effects [59] . Since the star tracker and gyro unit are located on the op- 
tical bench which is designed to hold the rigidity between instruments, the self-calibrating 
capability of gyro systematic errors by estimating gyro biases will be kept as long as the 
star tracker provides high frequency observations. The failure to achieve the star tracker 
measurements at the anticipated rate probably made the MSX attitude determination 
suffer from the MSX gyros systematic errors other than true biases. The GLAS requires 
the precise attitude/pointing determination in the high latitude regions such as Greenland 
and Antarctica (|<5| > 60°). Since the Sun and the Moon, the possible cause of the glint, 
are located near the ecliptic plane (|<5| « 23.5°), the effect of the Sun and the Moon light 
on the star tracker measurement should be minimal and uninterrupted tracker data rate 
is expected. Consequently, the Kalman filter that estimates both the attitude quaternions 
and the gyro bias simultaneously would provide the auto-calibration feature for the gyro’s 
systematic errors, if the HRG of the GLAS suffers from problems similar to the MSX gyro. 
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Chapter 7 


IMPLEMENTATION 

CONSIDERATIONS 

7.1 Standards 

To facilitate generation of the laser pointing direction, the attitude should be expressed 
using reference frames common to those used in description of the orbit position. As 
discussed by Rim and Schutz [35], the POD process will provide the position of the space- 
craft/instrument center of mass in the ICRF which is the specific implementation of the 
CRF. The determination of the transformation matrix between the spacecraft-fixed axes 
and the CRF depends, in part, on the adopted star catalog used with the analysis of star 
camera data. In mid 1997, the Hipparcos star catalog was released containing more than 
100,000 sources with position precision better than one mas (see Table 3.5). Although 
additional analysis will be required, it is expected that the Hipparcos catalog will be the 
foundation of the GLAS PAD.* Additional stars were available with the Hipparcos release, 
known as the Tycho catalog, but with precision at the 30 mas level. Nevertheless, the Ty- 
cho contribution contains about one million stars. The Hipparcos source coordinates will 
be defined for an epoch corresponding to J2000 and the reference system will be realized 
through the ICRF. The proper motions in the catalog are better than 1 mas/yr. The Hip- 
parcos catalog will link the optical sources (used by Hipparcos) to the radio sources (used 
by the IERS) at a level better than 1 mas [27]. Further discussion is given by McCarthy 
[ 31 1 - 

Conceptually, the determination of the laser spot location will be determined in the 
ICRF using the spacecraft center of mass ephemeris from the POD and the transforma- 
tion matrix from the PAD to determine the measured height vector. This height vector 
will be transformed into the ITRF using the IERS measured UT1, as well as precession 
and nutation parameters. The same transformation between the ICRF and the ITRF is 
required in the POD process. The transformation into the ITRF from the ICRF is well 
established at the milliarcsecond level. 

*The Hipparcos and Tycho catalogs provide high-quality astrometric and photometric data for virtually 
all SKY2000 stars in Version 2 and later versions. The SKY2000 Master Catalog (Version 2) was used in 
practice as the foundation of the GLAS PAD. 
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Table 7.1: Overall attitude determination accuracy goal 



roll(arcsec) 

pitch (arcsec) 

1. Attitude Determination 

0.47 

0.47 

2. Velocity Aberration 

0.03 

0.03 

3. Star Position Accuracy 

0.03 

0.03 

4. Ephemeris 

0.01 

0.01 

RSS(2,3,4) + 1 

0.52 

0.52 


7.2 Ancillary Inputs 

Ancillary inputs include a subset of a master star catalog/ and a luni-solar-planetary 
ephemerides, such as DE-200 or later [55]. The latter complements the star catalog by 
enabling identification of solar system bodies in camera images. For the generation of 
the transformation between the spacecraft-fixed axes and the ICRF, the velocity aberra- 
tion corrections of the observed stars are needed (see Section 3.3.3 and Table 7.1). The 
ephemeris of the GLAS satellite determined in the POD processing will be known with 
accuracy better than the real time flight operations ephemeris. In addition, to expedite 
the PAD process, access to the real time attitude determination from the flight operations 
segment is required. This information will reduce the search time to identify stars in the 
star camera images by providing a priori attitude information of modest accuracy. 

The POD process requires real time attitude information as well as a modest accuracy 
a priori ephemeris. The present algorithm development will pursue independent POD and 
PAD processes to simplify software validation, even though these processes (including the 
laser pointing process) share common aspects. The common parameters will be made 
available through shared files. 

7.3 Accuracy and Validation 

An estimate of the attitude determination accuracy is shown in Table 7.1 [60]. While 
the star catalog is expected to introduce no significant errors into the estimation process, 
the star camera will be the dominant source of error. The RSS of roll and pitch, which 
will affect the laser pointing direction, is 0.85 arcsecond. The overall RSS of all the error 
sources for the laser pointing determination is expected to be 1.20 arcsecond (la) (see 
Table 1.1). The simulation in this document did not include unmeasured errors, stemming 
from the unsensed bench vibration and the unsensed angular distortions. Even though 
the pointing errors shown in the plots in Chapter 6 will increase slightly when these 
unmeasured errors are taken into account, the resultant la value should be smaller than 
1.20 arcsecond considering the magnitude of the unmeasured errors and our simulation 
results. 

Detailed discussion for the attitude validation is contained in the Calibration/ Validation 
Plan [61]. In summary, a first order validation of attitude will be conducted by compar- 
isons with the real time determination. This approach, however, is limited because the 
real time requirements are about 10-20 arcseconds (la), compared to 1 arcsecond (la) 

^The master star catalog is SKY2000 catalog. The subset is GLAS PAD mission star catalog. 
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required for science analysis. As described in this plan, validation of the attitude at the 
1 arcsecond level is inherently linked to laser pointing. Experiments to validate the laser 
pointing will simultaneously assess the PAD. 

7.4 Computational: CPU, Memory and Disk Storage 

The requirements for the product of the PAD, the 3x3 transformation matrix between 
spacecraft-fixed axes and the ICRF, will require about 300 MB/day. Alternatively, about 
150 MB/day would be required if 4 quaternion elements were used at each time instead 
of 9 matrix elements. Preliminary tests show that an HP735-level machine can support 
the required attitude determination; however, further tests will be conducted to provide a 
more definitive estimate for CPU, memory and disk storage. 

The approximate CPU times for one orbital period of data processed with the EKF 
algorithms on a HP model 735/125 Unix machine were given in Table 5.2. From the 
table, the processing time for one day attitude of data is estimated to be approximately 
45 minutes for the EKF using the QUEST observation model and 90 minutes for the one 
using the UVF model. The batch algorithm is expected to require about 135 minutes for 
the same amount of data. The overall PAD process will require additional CPU time for 
the SRS data processing discussed in Chapter 6 and the data preprocessing. 

7.5 Sensor Failures 

Although the PAD study has been done with a star tracker in the instrument optical 
bench, there will exist two additional CCD star trackers (Ball CT-602) mounted on the 
spacecraft structure. The primary purpose of the spacecraft star trackers is to provide 
data for the real time attitude determination. The spacecraft star tracker data will be 
available on the ground and the additional star measurement data can be used for PAD 
in post-processing. Since the BDs of the CT-602 star trackers will not be parallel to that 
of the instrument star tracker, the contribution of this star tracker to the determination 
of the laser pointing direction is limited by the amount of the angle between star trackers’ 
BDs. However, in case the instrument star tracker malfunctions, the Ball CT-602 star 
trackers’ data must be used to produce the best laser pointing information. 

This idea can be extended when the optical bench star tracker loses data temporarily 
due to solar and lunar eclipsing with no permanent failure. Instead of totally depending 
on the CRS image stability as discussed in Section 6.1.5, it might be possible to use the 
attitudes obtained by the CT-602 star trackers whose BDs are 30° away from the optical 
bench star tracker. The calibrated data of the CT-602 trackers in terms of the instrument 
star tracker will provide 10 Hz attitude determination for the PAD with degraded accuracy. 

Meanwhile, the situation is not the same for gyro failure since there is no alternative 
system in the spacecraft. There are two possible ways to generate angular rates without 
the actual rate measurement even though the accuracy level is in doubt. The first choice 
is to use the Euler equation (Equation 2.18) to compute the spacecraft angular rates. For 
this case, the numerical or analytical expression for external torques must be modeled as 
a function of time as well as a function of the position and attitude of the spacecraft. 
The other choice is already mentioned in Section 5.4 when the EKF simulation results 
were discussed. The EKF without using gyro data showed equal or better accuracy in 
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attitude determination. However, it was indicated that the real spacecraft motion may 
contain high frequency noise which can be detected only by gyros. Even though the 
replaced system/method to the original PAD system will not determine the attitude to the 
required accuracy level, preparation of alternatives to produce the most accurate attitude 
determination as a backup plan is prudent. 
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Appendix A 


ICESAT MISSION OUTLINE 


The Ice, Cloud and land Elevation Satellite (ICESat) was launched on 13 January 2003. 
The Geoscience Laser Altimeter System (GLAS) instrument onboard ICESat made its 
first laser elevation measurement of the Earth on 21 February 2003 and its last on 11 
October 2009. The three lasers employed by GLAS did not perform as long as expected, 
and following the failure of Laser 1 on 5 March 2003 the ICESat mission was modified to 
meet the requirement for capturing a multi-year time series of ice sheet elevations [40]. 
For the modified mission scenario, the spacecraft entered a 91-day repeat science orbit 
(compared to a planned 183-day repeat) and the lasers were activated for about 33 days 
of this 91-day repeat, two or three times per year. This campaign mode operation is 
summarized in Table A.l, and other significant parameters and events are listed in Table 
A. 2. 


Table A.l: ICESat Laser Operation Campaigns 


Campaign 

Year 

Day of year 

Calendar Dates 

Number of 
days (d) 

Repeat 
orbit (d) 

Repeat tracks 1 

Lla 

2003 

051-088 

20 Feb-21 Mar 

37 

8 

001-072 to 
006-023 

L2a 

2003 

268-277/ 

277-322 

25 Sep-4 Oct/ 
4 Oct-21 Nov 

54 

8/ 

91 

028-088 to 029-100/ 
1098 to 0421 

L2b 

2004 

048-081 

17 Feb-21 Mar 

33 

91 

1284 to 0421 

L2c 

2004 

139-173 

18 May-21 Jun 

34 

91 

1283 to 0434 

L3a 

2004 

277-313 

3 Oct-8 Nov 

37 

91 

1273 to 0452 

L3b 

2005 

048-083 

17 Feb-24 Mar 

35 

91 

1258 to 0426 

L3c 

2005 

140-174 

20 May-23 Jun 

34 

91 

1275 to 0421 

L3d 

2005 

294-328 

21 Oct-24 Nov 

34 

91 

1282 to 0421 

L3e 

2006 

053-087 

22 Feb-28 Mar 

34 

91 

1283 to 0424 

L3f 

2006 

144-177 

24 May-26 Jun 

33 

91 

1283 to 0421 

L3g 

2006 

298-331 

25 Oct-27 Nov 

33 

91 

1283 to 0423 

L3h 

2007 

071-104 

12 Mar- 14 Apr 

33 

91 

1279 to 0426 

L3i 

2007 

275-309 

2 Oct-5 Nov 

34 

91 

1280 to 0421 

L3j 

2008 

048-081 

17 Feb-21 Mar 

33 

91 

1282 to 0422 

L3k 

2008 

278-293 

4 Oct- 19 Oct 

15 

91 

1283 to 0145 

L2d 

2008 

330-352 

25 Nov-17 Dec 

22 

91 

0096 to 0423 

L2e 

2009 

068-101 

9 Mar- 11 Apr 

33 

91 

1286 to 0424 

L2f 

2009 

273-284 

30 Sep-11 Oct 

11 

91 

1280 to 0084 


1 There are 119 tracks in the 8-day orbit and 1354 tracks in the 91 -day orbit. Cycle 
numbers are included for the 8-day repeat periods. 
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ICESat laser campaigns are designated by a laser number (LI, L2 or L3), followed by 
a letter in the sequence of operation. Following campaign L2f, attempts to restart any of 
the lasers were not successful. The spacecraft was put through a series of engineering tests 
in early 2010. De-orbit maneuvers were carried out in June and July 2010. The spacecraft 
was passivated on 14 August and reentered the Earth’s atmosphere on 30 August 2010 
over the Barents Sea northeast of Norway. 


Table A. 2: Significant ICESat Parameters and Events by Campaign 


Cam- 

paign 

Year 

Day 

of 

year 

S/C 

orient- 

ation 1 

Start 

Beta’ 

Angle 

o 

End 

Beta’ 

Angle 

(°) 

Start 

Laser 

Infrared 

Energy 

(mJ) 

End 

Laser 

Infrared 

Energy 

(mJ) 

Mean 
footprint 
major 
axis (m) 

Day of year — 
comments 

- 

2003 

013 

- 

- 

- 

- 

- 

- 

0 1 3 - launch 

Lla 

2003 

051-088 

-Y/+X 

-45 

-32 

72 

51 

149 

080 - yaw flip 
085 - safe hold, 
adjust temperature 

L2a 

2003 

268-277/ 

277-322 

+Y 

51 

69 

80 

55 

100 

277 — orbit change 
286 — laser 
temperature anomaly 
287, 302 - adjust 
temperature 
311 - GPS solar 
flare anomaly 

L2b 

2004 

048-081 

+Y 

54 

40 

57 

33 

90 


L2c 

2004 

139-173 

-X 

13 

-4 

33 

5 

88 

142-147 - adjust 
temperature 

L3a 

2004 

277-313 

-Y 

-48 

-58 

67 

62 

56 

293 - adjust 
temperature 

L3b 

2005 

048-083 

-Y 

-56 

-45 

68 

54 

80 

054 — suspected 
amplifier bar drop, 
begin footprint 
anomaly 2 
068 — suspected 
amplifier bar drop 

L3c 

2005 

140-174 

+X 

-20 

-4 

49 

44 

55 


L3d 

2005 

294-328 

+Y 

51 

63 

43 

39 

52 


L3e 

2006 

053-087 

+Y 

62 

48 

38 

30 

52 


L3f 

2006 

144-177 

-X 

20 

4 

30 

30 

51 

149 - Energy jump 
up 2m J 

L3g 

2006 

298-331 

-Y 

-44 

-54 

30 

24 

53 

310 -begin ITRF 
2005 

L3h 

2007 

071-104 

-Y 

-60 

-47 

24 

21 

56 


L3i 

2007 

275-309 

+Y 

32 

46 

22 

20 

57 


L3i 

2008 

048-081 

+Y 

74 

62 

20 

16 

59 


L3k 

2008 

278-293 

+X 

-28 

-32 

18 

12 

52 

289 - Energy drop 4 
mJ 

L2d 

2008 

330-352 

-Y 

-45 

-53 

8 

4 

- 

343-344 - adjust 
temperature, energy 
up 5 mJ 

L2e 

2009 

068-101 

-Y 

-71 

-59 

6 

2 

- 

094-095 - adjust 
temperature 

L2f 

2009 

273-284 

-X 

20 

25 

4 

2 

- 


- 

2010 

242 

- 

- 

- 

- 

- 

- 

242 — reentry 


1 The spacecraft is said to be in “Sailboat” mode for ±Y orientations and in “Airplane” 
mode for ±X orientations, where the direction indicates the solar panel orientation with 
respect to the spacecraft velocity using the GLAS coordinate frame. 

2 The footprint diameter during L3b changed from a mean of 54 m (day of year 048-053) 
to 84 m (055-068). The reason for the larger footprint size during the latter part of the 
campaign is unknown, although a suspected amplifier bar dropout occurs near the event. 


Appendix A was written by Tim Urban at Center for Space Research, The University of Texas at 
Austin. 




Appendix B 


PROCESSING METHODS FOR 
POINTING DETERMINATION 


The pointing determination algorithm described in Chapter 6 was developed based on 
the intended SRS behavior. In other words, the algorithm was a pre-launch processing 
design corresponding to the pre-launch SRS design. After launch, there were significant 
problems with the SRS system and the pointing determination algorithm was adapted to 
on-orbit SRS performance. For example, there were significant problems with the Laser 
Reference Sensor (LRS). It was deactivated after two weeks of intermittent operation [38]. 
Before the second laser campaign, in September 2003, an LRS software change resulted 
in continuous but degraded LRS operation. The degraded LRS measured CRS and laser 
spots continuously, but stars were only acquired and tracked during eclipse. No stars were 
measured in daylight. The LRS was designed to tie the LPA reference frame to the 1ST 
reference frame. Without continuous LRS star measurements, the LRS could not be used 
to transform laser measurements from the LPA reference frame into the 1ST reference 
frame. Significant modifications were made to the algorithm described in Chapter 6, 
particularly for the LRS. The modified algorithm was described in detail in reference [6]. 

B.l 1ST Bracket Motion 

Degraded LRS star measurements limited the SRS self-calibration capability described in 
Section 6.3.1. This place more emphasis on LRS measurements of the CRS and laser spots. 
These measurements included evidence of orbital variations of the spot positions. The 
CRS was attached to the 1ST mounting bracket. CRS spot motions reflected 1ST bracket 
motions. Smaller but similar motions of the laser spot suggested thermal variations of 
the GLAS optical bench due to sunlight. Early in the mission, the CRS and laser spot 
variations were incorporated into the processing algorithm [6]. 

The CRS signal was generated using green wavelength light from the GLAS laser. The 
energy of this light source decreased rapidly from January 2004. From the fall of 2004, 
the LRS was no longer able to measure the CRS spot. A batch-EKF correction method 
was developed to replace the CRS spot measurements [4], It compared a batch method 
attitude with an EKF attitude. The batch method attitude was estimated using gyro 
measurements alone, so was independent from the 1ST bracket motion. The EKF attitude 
was estimated using 1ST measurement so included the bracket motion. The difference 



between the two attitudes represented the 1ST bracket motion. 1ST measurements were 
corrected using the batch-EKF results. 

B.2 Time-tag Problems in 1ST and SIRU Data 

Several issues were found in the 1ST and SIRU time-tags: a 0.1 seconds 1ST time-tag 
shift, periodic SIRU rate spikes, SIRU data compression and bias. These problems caused 
anomalies in the attitude and pointing solutions. Analysis of the 1ST and SIRU data 
revealed the root causes and corrections were applied during reprocessing, as discussed in 
reference [2]. 

B.3 1ST Sun-Blinding 

ICESat’s two primary attitude modes were sailboat mode and airplane mode. Each pri- 
mary mode had two sub-modes separated by a 180° rotation. The primary mode was 
selected based on the beta prime angle (the angle between the orbit plane and the sun). 
Airplane mode was used when the absolute value of the beta prime angle was less than 
33°. Sailboat mode was used otherwise. The 1ST was periodically blinded by the sun 
during airplane mode campaigns. For example, the sun-blinding data gaps on May 20, 
2005 (beta prime angle —20°) were approximately 400 seconds long. On June 23, 2005 
(beta prime angle —4°) the data gaps were approximately 700 seconds long. The 1ST was 
also intermittently blinded by the moon, usually for less than 100 seconds. 

When the 1ST was blinded and star measurements were not available, attitude esti- 
mation was based on gyro data alone. The gyro data included errors from several sources 
including bias instability, scale factor error, and sense-axes misalignments. Attitude er- 
rors could grow to tens of arcseconds when using gyro data alone. There were several 
alternatives for improving attitude accuracy during 1ST blinding. One method was to 
use measurements from the BSTs, but this introduced uncertainty because the BSTs were 
pointed 30° away from the zenith. A method based on batch processing to estimate ad- 
ditional gyro parameters such as misalignments and scale factor errors resulted in better 
attitude estimates during data gaps [39]. A slight modification of the variable interval 
batch method [3] also provided good attitude estimates when the batch interval included 
approximately 100 seconds of star tracker data at the beginning and end of the data gap. 

During 1ST data gaps, batch processing was used instead of the EKF. For verification, 
the batch estimates were compared to the EKF estimates for the BSTs [39]. The batch 
estimates were also compared to EKF estimates before and after the data gaps. This was 
particularly useful when the beta prime angle was near 0° and both of the BSTs were 
blinded as well, for example during June 2004. 

B.4 Distortion Correction 

Star measurements were represented in the GLA04 data as unit vectors in a tracker coordi- 
nate frame. Two methods were used to express the unit vectors. BST1 and BST2 output 
two scaled angles kO ^ and kd v where the scaling factor k converted from radians to arcsec- 
onds. To convert the angles to a unit vector y, an intermediate vector x = [x\ X 2 x^} T was 
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calculated with X 3 defined equal to one, x\ defined as the horizontal coordinate h, and X 2 
defined as the vertical coordinate v, so that x = [h v 1} T where 


h = X1/X3 = x\/l = tan Oh 
v = X2/X3 = X2/I = tand v 

The unit vector y was then the normalized version of x 

y = [h v 1 } T /(h 2 + v 2 + 1 ) 1/2 (B.l) 

The native outputs of the 1ST and LRS were kh and kv. The 1ST and LRS performed 
the trigonometric calculations h = t&n 9h and v = tan#*, internally. Because 9 h and 9 V 
were small angles, kh and kv were nearly equal to k6h and k 9 v , but confusing the two 
types of output produced errors on the order of thirty arcseconds across an eight degree 
field of view. From GLA04 Release 24 (2005.12.21) onwards the 1ST values were converted 
from kh and kv to kQh and kQ v . If LRS star measurements are to be used, it is helpful 
to validate the LRS coordinates using a case where the 1ST and LRS are simultaneously 
measuring the same star. There is a ninety degree rotation between the 1ST and LRS 
coordinate frames. There is also a release dependent scaling factor for the LRS outputs, 
defined in the GLA04 release notes. 

Distortion corrections were applied to the measured star unit vectors before use. The 
corrections were estimated using a least squares fit of third-order polynomials to a large 
sample of measurements. An iterative method was used to accumulate the factors A = 
H T H and B = H T y for the least squares solution, where y = [y\ 2 / 2 ] T was an effective 
measurement. The iterations were initialized with Aq = O 20 X 20 and Bq = (Loxi- Given a 
measured unit vector u = [ui U 2 {u\ + U 2 ) 1 /' 2 ] 7 ' and a predicted unit vector v = [v\ V 2 (v 2 + 
v 2 ) 1 ^ 2 ] T expressed in the tracker frame, the measurement components were y\ = u\ — v\ 
and 2/2 = U 2 — V 2 - Two measurement models 2/1 = h\x + £\ and 2/2 = ^2 x + £2 were defined 
using a twenty-component vector of coefficients x and third-order polynomials 

hi = [1 U\ U2 u\ U\U2 U2 Ui U\U2 U1U2 «2 Olxio] 

9 909 OQ-i 

h2 = [Olxio 1 u \ U2 U1U2 U2 u\ U1U2 U1U2 U2 ] 

Aj and Bj were accumulated by 

Aj = Aj—i -{- hi hi T /^2 ^2 
Bj = Bj - 1 + h^yi + h^y2 

The estimated coefficients x were calculated using the least squares normal equation 

x = (iL T iL) - 1 Lf T 2/ = A~ X B 

The corrected unit vector w for the measured unit vector u was then 

r 9 o R o o q, r -, r T' 

W 1 = Ui — [1 Ul U2 U1U2 u 2 III Uiu 2 U1U2 U2\[Xi...Xiq\ 

r p 9*^9 9 o-, r -t r r 

W2 = U2 — [1 Ul U2 Ui UlU2 U 2 Ui Ui‘U 2 UlU 2 U 2 \ [ill • • -® 2 oJ 
w = [rci W2 (w\ + 

Distortion correction magnitudes were up to 2.5 arcseconds for the 1ST and up to 5 arc- 
seconds for the BSTs. Detailed descriptions of the distortion corrections are available in 
a reference [51]. 
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B.5 Bad Star Removal 


Predicted and measured star positions were compared and used to form star position 
residuals. The residuals were a direct measure of the attitude estimate accuracy. Daily 
mean values of the residuals were used to monitor estimation performance and in several 
cases were used to detect bad stars. The first case occurred during the fall 2003 laser 
campaign. The mean residual values were large on several days and were traced to a 
single star with large residuals. Excluding the star from the attitude estimation process 
resulted in normal mean residuals values. Here the definition of a bad star is a star with 
large residuals. 

Figure B.l shows the individual star position residuals over one orbit from 2008 day 
50. The residuals are cutoff above twenty arcseconds because all such measurements were 
automatically rejected in an early stage of the processing. The residuals were generally 
less than ten arcseconds, but strong peaks are apparent near 10,800 seconds and 11,450 
seconds. The individual stars that were being measured during the peaks were designated 
as possible bad stars or suspects. Means and standard deviations for the position residuals 
of individual suspects were then used to classify them as normal or bad. The bad stars were 
included in a list of stars that were automatically rejected in an early stage of processing. 
Figure B.2 corresponds to Figure B.l but show the results after bad stars were detected 
and rejected. 

A relatively small number of stars, less than one percent, were found to have statisti- 
cally significant position and brightness biases, primarily due to nearby bright stars. These 
stars are discussed in detail in references [17] [52] [53]. 
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Figure B.l: Star position residual in day 50 Figure B.2: Star position residual in day 50 
of year 2008 of year 2008 after bad star removal 
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Appendix C 


LPA IMAGE 
CHARACTERISTICS 


The final PAD product was referred to as ANC09. It was generated daily and sent to 
GSFC where it was combined with orbit and range information to calculate earth surface 
elevation. ANC09 contained a time-series of PAD pointing vectors that were estimated 
using the algorithm described in Appendix B. It also contained parameters describing the 
laser images measured by the LPA. The parameters included centroid position, orientation, 
major-axis, eccentricity, and total intensity. Orientation was represented by two angles: 
angle of major-axis from LPA x-axis, and angle of major-axis from topocentric north. The 
parameters were calculated using algorithms in Source Extractor, a program designed to 
catalog objects in astronomical images [9]. 

LPA data consisted of images with 80x80 pixels at 40 Hz. The Stellar Reference System 
was designed to capture a portion of the laser energy within a 20x20 portion of the LPA. 
The image pixel intensity values represented the energy distribution in the laser beam. 

LPA images of Laser 1 initially had an elliptical center with an irregular shape at one 
edge as shown in Figure C.I.* The irregular shape grew and detached from the central 
ellipse over about 30 days as shown in Figure C.2. Laser 2 laser images were elliptical like 
Laser 1, but without irregularities. An example of Laser 2 image is shown in Figure C.3. 
The ellipse orientation changes over time but the size and shape are almost consistent. 
Laser 3 laser images were roughly circular and constant over years of operation. Figure C.4 
shows an example of Laser 3 laser image. LRS image is similar to LPA image except for 
a 90° rotation due to their coordinate rotations to each other [5]. 

LPA characteristics were calculated for all 40 Hz data. For each laser image, 1/e 2 
rule was applied as a cut off value : pixels having intensity less than a cut-off value were 
ignored. Daily averages of LPA characteristics were calculated from these 40 Hz output. 
Then campaign averages were calculated from the daily averages in each campaign. Fig- 
ure C.5 shows changes of orientation (about LPA x-axis), major-axis and eccentricity of 
laser image, from LI (2003) to L3k (2008) campaigns. Figure C.6 has intensity informa- 
tion: total intensity, maximum intensity, and number of pixels whose intensity are greater 
than cut-off values. Figure C.7 shows centroid coordinates of LPA images: in x-axis, in 
y-axis, and the overall distance from the LPA origin, (0,0). As indicated in Table A. 2 
of Appendix A, the major-axis value changed abruptly on day 54 of 2005. It explains 

*LPA images are reconstructed using IDL from digitized numbers. 
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Figure C.l: Laser 1 LPA image (067/2003) 


Figure C.2: Laser 1 LPA image (082/2003) 



Figure C.3: Laser 2 LPA image (145/2004) Figure C.4: Laser 3 LPA image (283/2007) 
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Figure C.5: LPA image characteristics : Campaign averages 
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Figure C.6: LPA image characteristics : Campaign averages 
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Figure C.7: LPA image characteristics : Campaign averages 


abnormal jumps in major-axis and eccentricity plots for L3b campaign. The size of the 
image went to back down to previous range after L3b campaign. Daily variations of LPA 
characteristics in all laser campaigns are shown in Figure C.8, Figure C.9, and Figure C.10. 
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Figure C.8: LPA orientation, major-axis and eccentricity : daily averages 
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Figure C.9: LPA total intensity, maximum intensity, number of used pixels : daily averages 
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Figure C.10: LPA centroid location: in x-axis, in y-axis, from origin : daily averages 
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Appendix D 


PROPERTIES OF 
QUATERNION 


The four component quaternions are defined in Equation 2.1 with a constraint Equa- 
tion 2.2. They are often divided into the three component vector, q, and one scalar 
component, qq : 


9 = 


q 

94 


[91,92, 93, 94] T 


(P-1) 


The inverse quaternion is defined as : 



(D.2) 


This corresponds to the reverse rotation of the same amount of angle A 9 about the same 
rotation axis. 

Quaternions, as well as Euler angles, can be used to construct 3x3 rotation matrix. 
Unlike the Euler angles, which have 12 different sets for the rotation matrix, quaternions 
have a unique combination for the rotation matrix such as : 


C(q) 


9i - 92 - 93 + 94 2 (9 i 92 + 9394 ) 2(qiq 3 - q 2 94) 

2(9192 - 9394) — 9i + 92 - 9 3 + 94 2(q 2 q 3 + 9i94) 

2(9193 + 9294) 2(q 2 q 3 - 9194) ~qf ~ 92 + 9| + 94 


(P-3) 


It is easy to 
the rotation 


see that the change of signs for all quaternions simultaneously does not affect 
matrix. Inversely, 


9i = 

T (O 23 ~ C 32 ) 

4^4 


92 = 

4S (C31 “ Cl3) 


93 = 

4S (Cl2_C2l) 

(D-4) 


94 = ±-(Cn + C 22 + C 33 + 1) 1//2 , 
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where C\j is a matrix component in zth row and jth column. A computationally convenient 
algorithm has been developed as follows [45]. 


ql 

= \(l + 2C n 

- T ) 


ql 

= — (1 + 2C*22 

- T ) 


qi 

= -(I + 2C33 

~ T ) 

(D.5) 

ql 

= — (1 + 2C44 

-n 


CJ 

11 

bn 

= trace C = C n - 

+- C 22 + C 33 

(D.6) 


where 


and selects the q t which has the largest absolute value assuming it to be positive. The 
other three q/s can be obtained by dividing q t into the appropriate three of the following 
six equations: 


< 7 i <74 = ^ (C23 - C32) 

Q3Q4 = ^(Cl2 — C*2l) 

Q3Q1 = -^{Czi + Ciz) 


q 2 qi = -(C'31 - C13) 
<72 (ft = ^ (C23 + C32 ) 
Q1Q2 = |(Cl2 + C2I ) • 


(D.7) 


One important property is the simple form for combining sequential rotations, q' and 


q " , into an equivalent single rotation, q , by 


9 = 


1 

q'l 

-q'l 

1 — 

"5? 

-q'l 

q'l 

q'l 

q'l 

q'l 

-q'l 

q'l 

q'l 

l 

1 

-q'l 

-q'l 

1 


/ 

q ■ 


This is the quaternion composition that is usually written as the following form 

q = q <S> q 

using natural order, while some literatures prefer to write the same equation as 

/ „ // 
q = q © q 


(D.8) 

(D.9) 

(D.10) 


using historical order. Thus any set of g’s can be solved universally as a simple, nonsin- 
gular, bilinear combination of the other two. The matrix representation for quaternion 
composition is 

C(q) = C(q")C(q'). (D.ll) 

Therefore, many equations for quaternion computation are easily derived using the rotation 
matrix manipulations. 

From the definition in Equation 2.1, quaternions can also be parameterized in terms 
of the Euler angles : 

q(l,0) = [sin(0/2) 0 0 cos(0/2)] T 

q( 2,0) = [0 sin(0/2) 0 cos(0/2 )] T (D.12) 

q( 3,0) = [0 0 sin(0/2) cos(0/2)] T , 
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where 6 is a rotated angle and 1,2,3 are rotation axes employed sequentially. For example, 
the 3-1-3 Euler angles rotation can be written as 


9313(01,02,03) = 9(3,03) < 8 > q(l, 0 2 ) 0 9(3, 01 ) 
sin 62 cos(0i — 03) 
sin 02 sin (0i — 63 ) 
cos 02 sin(0i + 03) 
cos 02 cos (01 +03) 


(D.13) 


To convert quaternions to Euler angles, the rotation matrix must first be calculated 
from quaternions by Equation D.3. Then, the Euler angles are derived from the rotation 
matrix according to the chosen set of Euler angle sequences. 

To transform Euler angles to quaternions, the inverse process of the previous method 
can be used. A more convenient way for use in computer programming is as follows [22] : 


’ 94 " 


' 1 ' 

91 

= R 

0 

92 


0 

. 93 . 


_ 0 _ 


(D.14) 


where 


R = (cos ^0ii?4 + sin ^0ii? 3 )(cos \O 2 R 4 , + sin ^02-R2)(cos - 63 R 4 + sin (*-0 3 i?i) (D.15) 


and 


Ri = x 


^3 = X 
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